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Preface 


The  problem  of  navigation  through  fog  was  brought  tragically  to  the  fore  during 
the  summer  of  1999  when  John  F.  Kennedy  Jr.’s  plane  crashed  into  the  Atlantic  off 
the  coast  of  Martha’s  Vinyard.  The  key  factor  here  was  the  inability  to  determine 
the  vertical  direction  due  to  the  lack  of  visual  cues.  Less  catastrophic,  but  of  more 
common  concern,  are  frequent  delays  caused  by  fog  at  America’s  major  airports. 
It  is  estimated  that  one  major  air  package  carrier  alone  incurs  yearly  costs  mea¬ 
sured  in  millions  of  dollars  due  to  delayed  or  rerouted  flights  caused  by  airport 
fog  conditions.  Such  delays  would  not  be  necessary  if  alternatives  to  the  current 
methods  of  handling  fog  conditions  were  available. 

The  current  method  of  coping  appears  simply  to  involve  reducing  the  rate  of  air¬ 
craft  arrivals  and  departures  such  that  collisions  are  avoided.  However,  such  an 
approach  does  not  avoid  the  problems  encountered  by  John  Kennedy  Jr.,  nor  does 
it  curtail  the  possibility  of  collision  with  obstacles  on  the  runway.  (Incidences 
of  wildlife  appearing  on  runways  under  low  visibility  conditions  have  also  been 
noted  in  the  past.) 

Though  several  methods  for  improving  the  ability  to  see  through  fog  have  been 
attempted,  the  one  of  interest  in  this  study  involves  the  use  of  stereoscopic  vision 
devices  coupled  with  laser  illumination  and  realtime  image  deblurring  software. 
But  a  fundamental  question  remains  as  to  the  capabilities  of  this  technology  when 
confronting  real  world  fog  environments.  While  models  of  fog  scattering  prop¬ 
erties  exist,  the  propagation  models  for  use  of  these  models  have  often  relied  on 
Monte  Carlo  analysis  of  scattering  by  fog  particles.  This  analysis  is  usually  in¬ 
complete  —  it  cannot  adequately  address  the  spatial  and  angular  structure  of  the 
radiation  that  reaches  objects  within  the  fog  volume  and  the  return  signals  from 
objects  in  the  fog.  In  particular,  if  comer  reflectors  were  mounted  at  the  sides  of 
the  mnways,  could  such  reflective  objects  be  viewed  through  the  fog?  To  answer 
this  question,  the  current  analysis  provides  the  physical  model  of  one-way  prop¬ 
agation  through  the  fog  to  a  retroreflector  embedded  in  the  fog,  and  a  return  path 
back  to  camera  receivers  mounted  on  each  of  the  aircraft  wings. 
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Executive  Summary 


This  report  discusses  a  mathematical  analysis  for  assessing  the  capabilities  of  a 
system  to  navigate  through  fog.  This  study  was  conducted  as  part  of  an  FY99- 
funded  U.S.  Army  Research  Laboratory  (ARL)  Director’s  Research  Initiative.  The 
main  project  objective  was  to  determine  the  optimal  means  of  implementing  a 
system  designed  to  enhance  the  capability  of  seeing  objects  through  low  visibility 
atmospheres.  From  a  military  standpoint,  all  weather  forces  will  never  become  a 
reality  unless  aircraft  can  navigate  safely  to  landing  under  adverse  ground  visibil¬ 
ity  conditions.  Often  ground  conditions  are  the  limiting  concern  due  to  low-lying 
haze  or  fog.  In  the  majority  of  these  cases,  it  is  not  the  act  of  taking  off  that  is  of 
concern  but  that  of  landing  under  reduced  visibility.  Additionally,  from  a  civilian 
standpoint,  reduced  visibility  conditions  can  often  hinder  pilots  from  landing  at 
major  metropolitan  airports.  On  a  yearly  basis  millions  of  dollars  of  business  are 
lost  simply  due  to  aircraft  rerouting,  delays  in  landing  due  to  reduced  visibility 
conditions,  or  both. 

A  recent  patent  application  by  Wendell  Watkins  of  ARL  holds  the  potential  of 
solving  this  perennial  problem  through  the  use  of  active  stereo  illumination  and 
observation  of  a  landing  field.  This  system  envisions  use  of  a  near-infrared  laser 
mounted  on  an  aircraft  using  pulsed  beams  to  illuminate  reflectors  outlining  the 
runway.  A  portion  of  the  laser  light  is  presumed  to  penetrate  the  fog  and  be  re¬ 
flected  from  retroreflective  comer  cubes  located  along  the  sides  of  the  mnway. 
Cameras  mounted  on  each  wing  of  the  aircraft  are  then  used  to  acquire  an  image 
of  these  returning  signals.  The  received  signal  is  passed  through  deblurring  soft¬ 
ware  and  presented  to  the  pilot  in  a  headmounted  display  using  stereo  goggles. 

The  basic  physical  model  of  this  system  involves  propagating  a  Gaussian  laser 
beam  forward  through  a  fog  layer  to  a  retroreflector.  The  reflected  beam  then 
passes  through  the  fog  a  second  time  to  a  receiver  camera.  We  study  the  charac¬ 
teristics  of  the  return  path  beam  in  terms  of  positional  and  directional  variations. 
The  mathematical  simulation  of  the  propagation  process  utilizes  the  small  angle 
approximation  when  dealing  with  the  largely  forward-scattering  fog  aerosols  stud¬ 
ied.  This  small  angle  approximation  results  in  a  simpler  version  of  the  equation  of 
radiative  transfer  being  used  in  the  propagation  analysis. 

The  fog  aerosol  used  was  obtained  from  the  U.S.  Army’s  Electro-Optical  Systems 
Atmospheric  Effects  Library  Phase  Function  Database  (Tofsted  et  al.,  1997). 
Analysis  of  data  in  this  database  provided  a  summed  Gaussian  form  for  the  for¬ 
ward  and  backward  scattering  portions  of  the  phase  function.  The  phase  function 

iTofsted,  D.  H.,  B.  T.  Davis,  A.  E.  Wetmore,  J.  Fitzgerrel,  R.  C.  Shirkey,  and  R.  A.  Sutherland, 
EOSAEL92,  Aerosol  Phase  Function  Data  Base  PFNDAT,  ARL-TR-273-9,  U.S.  Army  Research 
Laboratory,  White  Sands  Missile  Range,  NM,  1997. 
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forms  were  then  then  coupled  to  a  Gaussian  beamwave  model  of  the  laser  source 
beam.  The  effects  of  the  retroreflector  material  were  modeled  as  a  spatial  clip¬ 
ping  coupled  with  a  nondiffraction-limited  spreading  (angular  convolution)  and 
reflection  of  the  incident  beam.  The  return  path  equations  were  then  assessed  us¬ 
ing  similar  equations  used  to  assess  the  outbound  beam.  The  combination  of  out¬ 
bound  propagation,  retroreflector  interaction,  and  return  beam  propagation  then 
allow  for  overall  analysis  of  the  propagation  problem. 

Approach 

To  model  this  propagation  problem  we  used  a  series  of  related  approximations 
as  well  as  processing  techniques  made  possible  by  the  approximations.  The  first 
step  involved  invocation  of  the  so-called  small  angle  approximation.  This  step 
simplifies  the  radiative  transfer  equation  (RTE)  and  permits  the  scattering  aerosol 
to  be  modeled  as  varying  only  along  the  2;  axis.  This  approximation  is  also  related 
to  the  use  of  a  sum  of  Gaussian  terms  expressing  the  angular  structure  of  the 
scattering  phase  function  for  the  aerosol  studied.  Using  the  small  angle  RTE,  a 
fourth-order  Fourier  transformation  in  the  off-axis  spatial  position  vector  [f  = 
{x,  y)]  and  a  two-dimensional  small-angle  vector  [a;]  is  performed  that  allows 
a  general  theory  (c.f.  Smirnov,  1964)^  to  be  used  in  solving  for  the  propagated 
radiance  as  a  function  of  position  and  direction  at  a  retroreflector  plane. 

This  solution  is  then  used  in  conjunction  with  an  equation  for  the  form  of  a  Gaus¬ 
sian  beamwave  in  this  four-dimensional  system.  The  net  results  of  this  analysis 
can  then  be  transformed  back  from  the  Fourier  domain  using  numerical  techniques 
to  study  typical  scenarios. 

Conclusion 

The  methods  described  here  are  also  useful  for  addressing  other  propagation  mod¬ 
eling  problems  and  questions  and  should  prove  useful  in  answering  questions 
about  the  capabilities  of  the  proposed  navigation  through  fog  device. 


^Smirnov,  V.  l.,A  Course  of  Higher  Mathematics,  Volume  II,  Pergamon  Press/Addison- Wesley, 
Reading,  MA,  1964. 


1.  Introduction 


The  performance  of  any  electro-optical  imaging  system  depends  on  the  medium 
through  which  signals  are  received.  As  optical  systems  are  continually  being  up¬ 
graded  for  special  purposes,  it  is  reasonable  to  augment  system  design  by  studying 
the  impact  of  environment  on  system  performance.  One  means  of  studying  this 
problem  is  to  simulate  the  system  its  environment  such  that  system-controlling 
parameters  can  be  modified  and  tested  against  a  variety  of  atmospheric  conditions 
one  can  reasonably  expect  to  encounter.  In  the  case  studied  here,  we  specifically 
want  to  look  at  optically  dense  fog  conditions  which  would  hinder  the  normal 
performance  of  most  sensor  systems.  This  study  was  conducted  as  part  of  an 
FY99-funded  U.S.  Army  Research  Laboratory  (ARL)  Director’s  Research  Ini¬ 
tiative  (DRI).  The  DRI  program  provides  a  means  of  funding  breakthrough  tech¬ 
nologies  at  a  startup  level  within  ARL. 

This  report  discusses  a  portion  of  this  overall  program  involving  a  mathematical 
analysis  of  the  capabilities  of  said  system.  The  main  project  objective  was  to  de¬ 
termine  the  optimal  means  of  implementing  such  a  system  for  landing  aircraft  (to 
include  helicopters)  under  low- visibility  ground  conditions.  From  a  military  stand¬ 
point,  all  weather  forces  will  never  become  a  reality  unless  aircraft  can  safely  land 
under  adverse  visibility  conditions  such  as  ground  fog  or  heavy  haze. 

To  facilitate  this  advanced  capability,  Wendell  Watkins  of  ARL  has  submitted  a 
patent  (Watkins  et  al.,  2000)  to  allow  vision  through  fog  via  the  use  of  active  stereo 
illumination  and  observation  of  a  landing  field.  This  system  envisions  use  of  a 
near-infrared  (IR)  laser  mounted  on  an  aircraft  fuselage  that  uses  pulsed  beams  to 
illuminate  reflectors  outlining  a  runway.  A  portion  of  the  laser  light  is  presumed  to 
penetrate  the  fog  and  be  reflected.  Cameras  mounted  on  each  wing  of  the  aircraft 
are  then  used  to  acquire  an  image  of  these  returning  signals.  The  received  signal  is 
passed  through  deblurring  software  and  presented  to  the  pilot  in  a  head-mounted 
display  using  stereo  goggles.  The  type  of  deblurring  kernel  used  and  the  effective¬ 
ness  of  the  deblurring  software  will  in  some  respects  depend  on  the  variability  of 
the  fog  aerosols  encountered,  but  may  be  anticipated  knowing  the  general  char¬ 
acteristics  of  the  fog  and  its  effects  on  propagation  of  the  Gaussian  laser  source 
beam  through  the  fog. 

In  this  report  we  study  the  characteristics  of  the  return  path  beam  in  terms  of  posi¬ 
tional  and  directional  variations.  The  mathematical  simulation  of  the  propagation 
process  utilizes  the  small  angle  approximation  since  the  fog  aerosols  studied  are 
largely  forward  scattering.  The  small  angle  approximation  means  that  the  majority 
of  the  energy  is  scattering  in  a  small  angular  cone  around  the  initial  propagation 
direction.  Use  of  this  approximation  results  in  a  simpler  version  of  the  equation  of 
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radiative  transfer  being  used  in  the  propagation  analysis.  The  fog  aerosol  used  was 
obtained  from  the  Army’s  Electro-Optical  Systems  Atmospheric  Effects  Library 
(EOSAEL)  (Shirkey  et  al.,  1987)  Phase  Function  Database  (PFNDAT)  (Tofsted 
et  al.,  1997).  Analysis  of  data  in  this  database  provided  a  sum  of  Gaussian  terms 
forming  the  forward-scattering  portion  of  the  phase  function. 

The  specific  propagation  problem  considered  here  is  the  degree  to  which  the 
backscattering  from  the  fog  can  be  mitigated  by  looking  at  the  return  signal  from 
a  source  laser  propagating  from  one  wing  of  the  aircraft  if  the  camera  viewing  the 
return  is  located  on  the  opposite  wing.  Such  a  separation  will  do  two  things:  (1) 
Use  of  stereo  goggles  will  allow  an  observer  to  better  distinguish  features  in  the 
vision  field;  (2)  By  locating  the  radiation  laser  source  on  the  one  wing  and  the 
receiver  on  the  opposite  wing  it  will  tend  to  reduce  the  amount  of  backscatter.  A 
third  advantage  results  from  image  processing  of  the  received  signals  from  both 
wing  cameras.  There  will  be  a  differential  change  in  the  amount  of  backscatter 
that  spans  both  received  images,  but  this  backscatter  signal  will  have  an  opposite 
gradient  on  each  (the  image  from  the  left  wing  will  see  more  backscatter  to  the 
right  side  of  the  image  and  visa-versa  for  the  opposite  wing  camera  images).  There 
could  be  advantages  in  differential  processing  of  the  images  obtained  from  each 
wing.  However,  to  determine  the  extent  of  these  effects,  one  would  need  to  model 
the  propagation  problem.  This  model  is  outlined  in  this  text. 


2.  Forward  Propagation  Through  Fog 


In  this  chapter  we  consider  the  problem  of  mathematically  modeling  the  single 
path  propagation  of  a  Gaussian  laser  source  beam  through  a  fog  aerosol  to  a 
retroreflector-type  object.  Following  the  solution  of  this  problem,  chapter  3  will 
deal  with  the  interaction  of  this  propagated  beam  with  a  retroreflective  surface 
embedded  in  the  fog.  Chapter  4  will  then  consider  the  propagation  of  the  reflected 
beam  back  through  the  fog  to  an  off-axis  camera.  Finally,  chapter  5  considers  di¬ 
rect  beam  energy  scattering  at  large  angles  -  the  diffuse  radiation  field.  This  latter 
calculation  is  critical  in  evaluating  the  relative  complicating  effects  of  the  diffusely 
scattered  energy  in  masking  objects  of  interest  in  the  field  of  view  within  the  fog. 

2.1  The  Radiative  Transfer  Equation 

To  model  the  propagation  of  laser  light  through  a  fog  layer  the  small  angle  ap¬ 
proximation  can  be  utilized.  This  technique  allows  for  the  solution  of  the  radiative 
transfer  equation  (RTF)  using  Fourier  methods,  permitting  the  calculation  of  both 
the  directional  and  positional  structure  of  the  forward  propagating  radiance  field 
at  a  plane.  This  approximation  assumes  the  majority  of  the  energy  involved  in  the 
radiative  transfer  process  is  only  flowing  into  the  forward  hemisphere. 

The  small  angle  approximation  involves  a  modification  to  the  radiative  transfer 
equation  from  a  fully  three-dimensional  form  to  a  single  hemisphere  form.  But  to 
discuss  this  change  we  must  first  consider  the  fully  three-dimensional  form  of  the 
radiative  transfer  equation,  given  by, 

n-'VI{s,h)+a  =  as  f  I{s,Cl')  P{n,Cl')  dCl' +{a-as)b{T{s),  A).  (1) 

J47r 

The  primary  quantity  generated  by  this  equation  is  a  radiance  field,  I{s,  Cl).  This 
field  varies  as  a  function  of  position,  denoted  by  the  three-dimensional  vector  s, 
and  direction,  denoted  by  the  unit  vector  Cl.  Cl  is  composed  of  directional  cosines 
oriented  along  the  x,  y,  and  2;  coordinate  axes. 

Radiance  is  measured  in  units  of  W/fm^-sr),  which  specifies  the  rate  of  energy 
passage  through  a  differentially  small  surface  into  a  differentially  small  solid  an¬ 
gle.  In  this  analysis  the  spectral  dependence  of  the  radiation  is  suppressed,  but  in 
general,  a  series  of  monochromatic  radiative  transfer  calculations  must  be  spec¬ 
trally  integrated  to  obtain  results  for  band-averaging  sensors.  For  laser  radiation, 
the  waveband  over  which  the  laser  operates  is  normally  so  narrow,  and  the  spec¬ 
trally  varying  characteristics  of  aerosols  are  normally  so  slowly  varying  with 
wavelength,  that  the  radiative  transfer  of  such  energy  can  be  considered  relatively 
monochromatic. 
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2.1.1  Terms  within  the  Radiative  Transfer  Equation 

The  first  term  on  the  left  hand  side  (LHS)  of  the  radiative  transfer  equation  (RTE) 
represents  the  rate  of  change  in  radiance  due  to  propagation  in  direction  Cl.  This 
term  is  thus  a  directional  derivative,  dl/ds,  where  ds  is  the  magnitude  of  ds,  a 
differential  change  in  position  in  the  direction 

The  second  term  on  the  LHS  represents  the  rate  of  extinction  of  the  radiance  field 
due  to  both  absorption  and  scattering,  where  (x  is  called  the  extinction  coefficient 
and  is  given  in  units  of  inverse  distance.  Relating  this  term  with  the  first  term 
on  the  LHS,  the  product  ads  is  thus  the  probability  of  interaction  between  the 
radiance  and  particles  embedded  in  the  medium  within  the  incremental  path  length 
ds.  For  fog  aerosols,  a  ranges  between  5  km“^  and  30  km~\  implying  several 
collisions  (virtually  all  involving  scattering  phenomena)  per  photon  per  kilometer 
travelled. 

Due  to  radiation  scattering,  energy  lost  from  the  radiance  field  due  to  extinction 
will  most  likely  reappear  as  scattered  energy  via  the  first  term  on  the  right  hand 
side  (RHS)  of  the  RTE.  This  returned  energy  re-enters  the  propagating  radiance 
field  by  scattering  from  different  directions  into  the  path  of  the  original  photon. 
The  nature  of  this  modified  directional  field  is  quantified  in  terms  of  the  scatter¬ 
ing  phase  function  (P(0,  D')),  which  describes  the  probabilistic  distribution  of 
scattering  based  on  both  the  incident  (Cl')  and  scattered  (Cl)  directions. 

The  rate  of  scattering  is  proportional  to  the  scattering  coefficient,  as,  given  once 
again  in  inverse  distance  units.  <Ts  is  related  to  the  extinction  coefficient  a  through 
a  quantity  called  the  single  scattering  albedo,  w,  which  quantifies  the  probability 
of  a  scattering  event  (as  opposed  to  an  absorption  event)  given  that  an  interaction 
between  a  photon  and  a  particle  occurs.  Hence,  we  normally  write, 

as  = ‘CO  a.  (2) 

The  final  term  on  the  RHS  of  the  RTE  represents  emission  of  thermal  radiation, 
which  depends  on  the  radiation  wavelength  (A),  the  temperature  of  the  atmosphere 
at  that  point  (T'(s)),  and  the  absorption  coefficient,  given  as, 

cTa  =  (1  -  'Co)  a.  (3) 

The  function  b{T,  A)  represents  the  radiation  emitted  by  a  perfect  blackbody  (im¬ 
plying  an  assumption  of  thermodynamic  equilibrium).  The  absorption  coefficient 
modifies  the  blackbody  emission  to  account  for  variations  in  emissivity  associ¬ 
ated  with  the  properties  of  the  aerosols  and  gaseous  constituents  of  the  medium. 
However,  we  will  ignore  emission  effects  in  our  analysis  becauseemission  influ¬ 
ences  are  insignificant  at  near-IR  wavelengths  for  normal  ranges  of  atmospheric 
temperatures. 

2.1.2  Solid  Angle 

The  integration  contained  within  the  first  term  on  the  RHS  of  the  RTE  is  re¬ 
quired  to  evaluate  the  energy  scattered  into  the  output  direction  of  interest.  Cl. 
This  integration  is  performed  with  respect  to  a  variable  known  as  solid  angle. 
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Solid  angle  is  used  to  evaluate  how  radiant  energy  spreads  out  following  passage 
through  a  point,  or  conversely,  how  radiant  energy  converges  toward  a  point  as  it 
propagates.  The  meaning  of  solid  angle,  and  similarly,  how  one  measures  solid 
angle,  can  be  thought  of  as  an  extension  of  the  means  for  defining  and  measuring 
two-dimensional  angles  in  Euclidean  space:  A  two-dimensional  angle  is  measured 
by  drawing  two  straight  lines  that  originate  at  a  central  point  and  extend  to  two 
different  points  on  a  circle  a  radial  distance  r  away.  The  angle  subtended  between 
these  two  lines  is  computed  by  dividing  the  length  of  the  arc  (a)  between  these  two 
points  by  the  radial  distance  (r):  a/r.  Solid  angle  is  a  three-dimensional  analog 
of  this  concept:  Instead  of  extending  two  lines  from  a  point  to  a  circle,  a  conical 
surface  is  extended  from  a  point  to  a  circular  region  on  a  sphere.  We  then  measure 
the  surface  area  within  this  region  on  the  surface  of  the  sphere  and  divide  by  the 
square  of  the  radius  of  the  sphere.  In  a  more  general  sense,  we  are  not  restricted  to 
a  circular  region  of  the  sphere,  but  rather  may  define  some  arbitrarily  connected 
area  on  the  surface  of  a  sphere  that  subtends  the  solid  angle.  The  boundary  of  this 
surface  can  be  connected  by  an  infinite  set  of  straight  lines  to  the  center  of  the 
sphere.  In  either  case  the  solid  angle  (Q)  associated  with  this  surface  is  evaluated 
by  dividing  the  surface  area  (A)  of  the  region  by  the  square  of  the  radius  (r^)  of 
the  sphere:  Q,  =  A/r^.  The  units  of  solid  angle  are  called  steradians  (sr),  though 
from  a  unit  analysis  standpoint  the  length  dimensions  appear  to  cancel. 

Because  of  the  definition  of  solid  angle  as  the  ratio  of  an  area  on  a  sphere  to  the 
radius  of  the  sphere  squared,  one  often  encounters  the  concept  of  a  “solid  angle 
sphere.”  This  term  refers  to  the  solid  angle  when  the  area  A  consists  of  the  entire 
surface  area  of  the  sphere  of  radius  r.  For  this  particular  case  we  have  A©  =  47rr^ 
and  thus 

n©  =  A©/r^  =  47r  sr,  (4) 

is  the  solid  angle  associated  with  a  sphere,  and  can  be  referred  to  as  a  solid  angle 
sphere.  At  the  other  extreme,  for  a  differentially  small  region  of  solid  angle  (dfi), 
there  will  be  essentially  only  a  single  propagation  direction  (17)  associated  with  all 
the  radiation  propagating  away  from  a  single  point  (s)  that  falls  within  that  solid 
angular  region. 

2.2  The  Small  Angle  Approximation 

The  small  angle  approximation  recasts  the  three-dimensional  RTE  into  a  form 
that  emphasizes  propagation  within  a  small  solid  angle  cone  about  the  forward 
direction.  However,  though  we  are  really  interested  only  in  propagation  near  the 
forward  direction,  we  shall  see  that  we  may  apply  the  results  of  this  analysis  to 
propagation  in  both  the  forward  and  backward  hemispheres,  but  the  results  will  be 
most  accurate  in  the  forward  direction. 

The  rationale  behind  the  small  angle  approximation  is  that  for  many  aerosols,  fogs 
in  particular,  the  majority  of  scattering  of  radiation  occurs  in  nearly  the  forward 
direction.  Coincidently,  this  behavior  is  termed  forward  scattering  (van  de  Hulst, 
1980b).  We  take  the  forward  direction  along  the  positive  z  axis. 
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In  general,  the  propagation  unit  vector  Cl  is  defined  in  terms  of  angular  variables 
0  and  <f),  as  illustrated  in  figure  1.  In  this  system,  0  is  a  polar/zenith  angle  and  <j)  is 
an  azimuthal  angle,  such  that. 

Cl  =  f2a;Z  +  +  0,zk 

=  cos^sin9i  +  sm.<f>sm9j  +  cos0k.  (5) 


Figure  1.  Angular  coordinate  system  using  two  angular  measures,  9  and  (p. 


(6) 


In  the  forward  hemisphere  the  triplet  (fia,  fly,  flz)  can  be  expressed  as, 

(fi.,  S2,.  n.)  =  — ) , 

where  -oo<a;x<oo,  -oo<a;y<oo  represent  components  of  the  vector, 

u;  =  oil  ?  +  Wy  (7) 


lying  in  the  (wx,  ojy)  plane.  For  propagation  into  the  backward  hemisphere  the 
sign  on  fl^  is  reversed. 

Examining  the  meanings  of  uJx  and  Uy,  we  consider  figure  2  which  shows  the 
relationship  between  9  and  (jj  =  \Q\  =  +  (*Jy- 

An  analysis  of  this  figure’s  geometry  reveals  that  for  differential  d9,  triangles  abc 
and  cde  are  similar,  leading  to  the  relations, 


cos{9)  = 


V^l  + 


dlL  ’ 


sin(0)  = 


u 


VT+a^ 


(8) 
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Figure  2.  Geometric  relationship  between  angle  6  and  tangent  U). 


Hence  the  integral  over  scattering  angle  in  the  RTE  can  be  transformed  into  an 
integration  over  and  w'  in  the  u  forward  scattering  plane: 


dn'  =  sine' de' d(f>' = 


u'dw'd(f)' 

(l+a;'2f/2 


(9) 


We  now  consider  the  way  in  which  the  phase  function  representation  is  also  af¬ 
fected  by  the  small  angle  approximation  and  the  use  of  u.  The  issue  of  back  scat¬ 
tering  will  also  be  addressed,  but  the  results  of  this  discussion  will  not  be  applied 
until  chapter  5  which  discusses  diffuse  radiation  propagation. 

For  the  moment  it  is  instructive  to  mention  that  purely  forward  scattering  may 
be  sufficient  to  handle  our  propagation  considerations;  Through  proper  filtering  it 
may  be  unnecessary  to  consider  diffuse  scattering  at  all.  In  particular,  techniques 
such  as  (1)  time  gating,  (2)  doppler  filtering,  and/or  (3)  spatial  filtering  may  elim¬ 
inate  large  scattering  angle  effects.  Assuming  one  or  more  of  these  is  used,  only 
a  small  region  in  the  u  plane  need  be  considered  as  containing  relevant  infor¬ 
mation.  Another  consideration  is  that  the  forward  scattering  portion  of  the  phase 
function  is  usually  very  large,  and  is  even  larger  with  aerosols  such  as  fog,  where 
large  mean  droplet  sizes  result  in  substantially  increased  probability  of  forward 
scattering. 


2.3  Phase  Function  Representation 

The  scattering  properties  of  the  aerosols  considered  in  addressing  the  problem  of 
navigation  through  fog  are  detailed  in  the  Army’s  Electro-Optical  Systems  Atmo¬ 
spheric  Effects  Library  (EOSAEL)  (Shirkey  et  al.,  1987)  Phase  Function  Database 
(PFNDAT)  (Tofsted  et  al.,  1997).  These  scattering  properties  were  evaluated  based 
on  particle  size  distributions  appropriate  for  fog  aerosols  and  a  Mie  scattering  cal¬ 
culation  needed  to  determine  the  directional  scattering  properties  of  the  aerosols. 
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These  scattering  probabilities  are  contained  in  a  phase  function.  The  fog  aerosols 
considered  were  developed  using  the  Deirmendjian  C.  1  aerosol  model  (Deirmend- 
jian,  1969),  as  applied  in  Shettle  andFenn  (1979)  to  fog  aerosols.  This  particle  size 
distribution  information  was  then  combined  with  complex  index  of  refraction  in¬ 
formation  for  liquid  water  and  used  in  A.  Miller’s  (1983)  AGAUS  Mie  scattering 
code.  The  resulting  phase  function  information  could  then  be  analyzed  and  inter¬ 
preted  in  its  forward  scattering  form  as  detailed  in  the  remainder  of  this  section. 

As  seen  from  equation  (1),  the  phase  function,  P{Q,  Q'),  determines  the  structure 
of  the  scattered  radiance  field.  This  function  represents  a  directional  probability 
distribution  of  scattering  in  any  particular  outbound  direction  (fi)  given  a  value  of 
incident  radiance  arriving  from  direction  Cl'.  The  normalization  condition  used  is 
that  when  integrating  a  phase  function  for  a  fixed  direction  of  incidence  (Cl')  over 
all  outbound  directions  (0©  =  47r  sr)  the  result  is  unity: 


1. 


(10) 


Another  property  common  to  all  phase  functions  is  based  on  the  reciprocity  princi¬ 
ple  (van  de  Hulst,  1980a):  The  phase  function  is  constrained  such  that  P{Cl,  fi')  = 
P{Cl',  Cl).  In  other  words,  all  phase  functions  must  be  symmetric  between  input 
and  output  directions. 

For  some  species  of  scatterers,  the  phase  function  may  depend  on  complete  knowl¬ 
edge  of  both  the  incident  and  scattered  propagation  directions  (e.g.,  snowflakes, 
raindrops).  But  for  hydrosols,  such  as  fogs,  the  particle  shapes  are  nearly  spher¬ 
ical  and,  hence,  symmetry  arguments  reduce  the  phase  function  dependence  to 
only  the  scattering  angle  Og.  Sometimes  this  dependence  is  written  in  terms  of 
—  cos  Og  ™  Cl '  ^ . 

For  the  particular  problem  of  interest,  scattering  through  heavy  haze/fog,  while 
models  of  the  phase  function  in  its  complete  form  exist  (Tofsted  et  al.,  1998),  a 
more  usable  representation  is  desired.  One  such  form  is  a  purely  Gaussian  phase 
function  (c.f.,  Tam  and  Zardecki,  1980;  Zardecki  and  Gerstl,  1986).  In  this  form 
the  phase  function  is  modeled  as  purely  forward  scattering,  using  a  Gaussian  func¬ 
tion,  modeled  here  as 

Pg{Aoj)  =  exp  tt  ,  (1 1) 


where 


\uj  —  u'l 


6g  w  Auj  =  |Au; 
and  a  is  an  inverse  width  parameter. 

A  useful  feature  of  the  Gaussian  phase  function  is  that  it  is  separable: 


(12) 


cP'  exp  [^— TT  =  oP  exp  j^— Aw^j  exp 


-CPAOJy 


(13) 


The  advantage  of  this  feature  is  seen  more  clearly  when  the  RTE  is  modified  for 
the  small  angle  approximation  and  transformed  into  the  Fourier  domain.  Also  note 


(14) 


that  Pg  exhibits  the  following  normalization  condition: 

OO 

27r  J  Pg{oJ,  a)u)duj  =  1, 

0 

such  that  the  integral  of  Pg  over  the  u  plane  normalizes  to  unity. 

However,  it  is  difficult  to  use  a  purely  Gaussian  phase  function  representation 
to  model  an  actual  phase  function  because  there  is  usually  more  than  one  criti¬ 
cal  scale  size  responsible  for  the  behavior  of  the  phase  function.  Typically,  phase 
function  profiles  exhibit  a  highly  peaked  and  narrow  forward-scattering  region, 
due  to  diffraction,  and  a  broader  secondary  peak  extending  to  50-60°.  These  fea¬ 
tures  are  illustrated  in  figure  3,  which  shows  plots  of  the  1.5  /^m  and  1.54  /xm 
wavelength  phase  functions  for  advection  fog. 


Scattering  Angle  0s 

Figure  3.  Advection  fog  phase  functions  for  1.5  fJim  and  1.54  flm  wavelengths. 

To  accurately  portray  the  phase  function  structure  it  will  be  necessary  to  formulate 
the  phase  function  in  a  form  which  is  both  accurate  and  amenable  to  computation. 
Both  goals  are  achievable  if  we  can  express  the  phase  function  as  a  sum  of  Gaus¬ 
sian  terms  for  each  hemisphere  individually  using. 


M 


Pf{Au;)  =  E  AnPgiAuj,  a^y, 

m=l 

(15) 

N 

Pti^LO)  =  Yl^ri  ^n)- 

(16) 

n=l 
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As  a  normalization  condition  we  require, 


M  N 

l='£Am+YlBn,  (17) 

m=l  n=l 


such  that  the  combined  Gaussian  components  from  both  hemispheres  integrate  to 
unity,  which  reflects  our  phase  function  normalization  condition. 

However,  there  are  two  complications  which  tend  to  cancel  each  other:  (1)  In 
assessing  the  fit  coefficients  for  the  Gaussian  series  in  both  hemispheres  we  note 
that  the  phase  function  approaches  a  nonzero  value  at  Os  =  90°;  (2)  The  phase 
function  integral  contains  the  additional  factor  (1  +  in  the  integration 

kernel: 


+ 


II 

rm(|a+J"|) 
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I,  (s,  u')  da' 

h  (s,  Cj")  du". 


(18) 


Here,  u"  denotes  integration  over  the  backward  hemisphere  where  the  phase  func¬ 
tion  peak  occurs  where  u"  =  —Co. 

In  this  form  the  factors  (1  +  a;^)"^/^,  appearing  in  the  denominators,  reflect  a  de¬ 
creasing  weight  with  angular  separation.  But  they  would  significantly  complicate 
the  problem  when  the  Fourier  transform  is  taken  (section  2.5).  To  avoid  this  situ¬ 
ation,  since  |a5|  is  small  due  to  the  small  angle  approximation,  we  may  substitute 
the  integration  with  respect  to  Co'  and  Co"  with  integrations  with  respect  to  Au)  and 
Ato  such  that 


^  Tr  Pf  (\^' )  _ 

I{s,  Cl')  P{Cl,  ClI)  dQ'  =  jj  3-7172  If 


-oo  ( 1  +  j 


+ 


oo 

//. 


l-b  Aa;' 


«2\ 


,  (19) 


) 


allowing  us  to  define  the  modified  phase  functions. 


PfiAJ)  = 


(H-Aa;'2)3/2> 


PbiAio")  = 


PfiAto") 

(H-Aa;"2)3/2’ 


(20) 


and  attempt  to  fit  Gaussian  series  to  these  functions.  Gaussian  approximations  to 
the  functions  Pf  and  Pb  are  shown  in  figure  4  for  the  forward  hemisphere  and  in 
figure  5  for  the  backward  hemisphere. 

In  the  forward  hemisphere  approximation  we  used  nine  Gaussian  components. 
In  the  backward  hemisphere,  we  used  eight  Gaussian  components.  Both  of  these 


Modified  Phase  Function  ^  Modified  Phase  Function 


4.  Forward  hemisphere  portion  of  phase  function  for  advection  fog  at  1.5  /Xm  wavelength. 


Figure  5.  Backward  hemisphere  portion  of  phase  function  for  advection  fog  at  1.5  (Jim  wavelength. 
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functions  appropriately  approach  zero  at  infinity,  but  the  series  fit  is  only  designed 
to  be  valid  for  scattering  angles  ranging  from  0—78°.  The  modified  phase  function 
approach  is  completely  accurate  in  the  forward  and  backscattering  directions,  but 
exhibits  increasing  error  for  off-axis  directions.  Since  we  are  primarily  interested 
only  in  propagation  along  the  main  axis,  this  should  not  be  a  critical  problem.  The 
Gaussian  expansion  coefficients  used  are  listed  in  table  1. 

Table  1.  Gaussian  expansion  coefficients  for  forward  and  backward  hemisphere  phase  function  representa- 
tions  of  the  1.5  /im  advection  fog  phase  function. 


m,  n 

Bn 

1 

0.3423 

29.38 

0.00009 

16.22 

2 

0.1144 

9.911 

0.0158 

1.7150 

3 

0.2676 

1.584 

-0.0463 

1.2810 

4 

0.1884 

0.7401 

0.0554 

0.8540 

5 

0.0377 

0.3741 

0.0051 

0.4076 

6 

0.0082 

0.2079 

0.0025 

0.1775 

7 

0.0034 

0.1243 

0.0016 

0.0800 

8 

0.0017 

0.0667 

0.00093 

0.0290 

9 

0.00107 

0.0290 

— 

— 

2.4  SmaU  Angle  RTE 

The  key  change  in  the  RTE  occurs  when  the  directional  variables  in  the  RTE  are 
replaced  by  functions  of  Q  and  Q' .  When  these  functions  are  substituted  into  the 
RTE,  it  may  be  written  as, 

OO 

— — H —  ^  -  —+a  I  =  vja  [[Pf  (IAo'  1)  f  dAw  .  (21) 

Note  that  ul'  =  aJ  -  such  that  the  phase  function  integral  appears  as  a  convo¬ 
lution. 

This  equation  can  be  further  simplified  by  assuming  that  in  general  Cj  will  be  small 
in  magnitude,  allowing  us  to  ignore  the  Vl  +  factors  appearing  on  the  LHS. 
This  simplification  leads  to  the  equation, 

OO 

—  OO 

2.5  Fourier  Transform  Operation 

The  key  advantage  of  writing  the  RTE  in  the  small  angle  form  is  that  the  angular 
integration  on  the  RHS  of  the  modified  RTE  is  in  the  form  of  a  two-dimensional 
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convolution  integral.  Convolutions  Fourier  transform  into  products  of  the  trans¬ 
forms  of  the  two  convolved  functions.  It  would  therefore  seem  expedient  to  Fourier 
transform  the  RTF.  The  transform  process  is  conducted  in  four  dimensions:  x,  y, 
tVx,  and  Uy.  Let  it  be  the  frequency  space  representation  of  f  =  xz  -f  yj,  and  let  if 
be  the  frequency  representation  of  u. 

To  perform  these  transformations,  Gaskill’s  (1978)  version  of  the  Fourier  trans¬ 
form  will  be  used.  In  this  version, 

OO 

F{i^x)  =  j  f{x)  exp  [-j2'KKa:x]  dx.  (23) 

—  OO 


OO 


f{x)  =  J  F{k,x)  exp[j2TTKxX]  dKx, 

— 

(24) 

with  j  =  F{kx)  is  the  Fourier  transform  of  the  real  valued  function  f{x). 

Gaskill’s  corresponding  Hankel  transforms  are  given  by 

OO 

G{k)  =  27r  g{r)  Jq  [2'KKr]  r  dr, 

0 

(25) 

OO 

g{r)  =  277^  G{k)  Jq  [27r«;r]  k.  dK, 

(26) 

0 


where  G{k,)  is  a  radially  symmetric  frequency  representation  of  the  radially  sym¬ 
metric  function  g{r). 

Given  these  forms  for  the  Fourier  and  Hankel  transforms,  we  consider  some  gen¬ 
eral  transform  rules:  First,  if  g{r)  has  the  Hankel  transform  given  by  G{k),  then 
g{r/h)  transforms  as 

H{!!(9}  =  |i'PGW.  (27) 

This  rule  is  useful  in  transforming  the  phase  function,  since 

H  |exp  1^— TT  I  =  exp  [— tt  77^]  .  (28) 

Hence,  we  can  write  the  Fourier  transform  of  the  phase  function  as, 


M 

Pfiji)  =  'P-{PM}  =  IZ  exp  [-7r77Va^  • 


m=l 


(29) 


Several  other  rules  relating  to  Fourier  transforms  are  also  of  significance:  If  / {x) 
has  the  Fourier  transform  F{kx),  then, 


T {f{x  di  xo)}  =  F{kx)  exp  [±j2TvxoKx] ;  (30) 

F{f'^Hx)}  =  {j2TrKx)'^FM,  (31) 
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where  represents  the  fcth  derivative  of  f(x)  with  respect  to  x.  Due  to  the 

similarity  between  the  Fourier  and  inverse  Fourier  transforms, 

T  {f{x)  exp  [±j27rxo«x]}  =  F{kx  T  ^^o);  (32) 

:P{(-j27ra)/(i)}  =  F<‘>(K,);  (33) 

A  final  rule  applies  to  convolutions:  If  f(x)  and  g{x)  Fourier  transform  to  F(k.) 
and  G{kx)7  respectively,  then  their  convolution, 

OO 

f{x)*g{x)=  j  f{P)9{x  -  P)dP,  (34) 

— OO 

transforms  as, 

r{f{x)*g{x)]  =  F(K.)GM.  (35) 

These  rules  are  used  in  transforming  the  small  angle  form  of  the  RTE.  Using  the 
derivative  rule,  we  transform  the  term  Ux  dl /dx  with  respect  to  x  as, 

Fx  3‘2‘Tri^x^xix,  (36) 

where  4  is  the  transform  of  I  with  respect  to  x.  The  resulting  j27r  factor  is  used 
in  transforming  with  respect  to  Ux, 

=  -c.  (37) 


where  Ixu^  is  th®  transform  of  /  with  respect  to  both  x  and  ujx.  Completing  this 
derivation,  we  find, 


—  —Kr 


drfx 


=  — 


dl 

drix 


(38) 


where  /  is  the  fourth  order  Fourier  transform  of  the  radiance  function.  A  similar 
procedure  is  followed  to  obtain  the  transform  of  the  term  Uy  dl/dy. 

The  convolution  rule  is  used  when  transforming  according  to  w  on  the  scattering 
integral  term.  In  addition,  we  use  linear  properties  of  the  operators  to  reverse  the 
order  of  operation  similar  to  the  procedure  used  above  where  the  operator 
was  switched  with  the  differential  operator.  Here,  following  the  operation, 

the  transformed  phase  function  is  factored  out  of  the  transformation  step  since 
it  is  considered  a  constant  with  respect  to  the  x  and  y  dimensions: 


PfilLO-  u'\)  I  (f*,  w',  z)du'\=  Pf  /. 


(39) 
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The  net  result  of  the  Fourier  transform  process  is  that  the  RTF  can  now  be  ex¬ 
pressed  as. 


d  d 


waPf 


1  =  0, 


where  the  dot  product  on  the  LHS  is  an  abbreviation  for, 


(40) 


^  d  d  d 

07]  orjx  dr]y 


(41) 


The  notation  in  equation  (40)  uses  a  singly  carated  variable  to  symbolize  a  dou¬ 
bly  transformed  function  (in  particular,  and  Uy  transformed),  while  a  doubly 
carated  variable  has  been  quadruply  transformed  (x,  y,  and  Uy).  Solving  for 
the  radiance  in  Fourier  space  is  viewed  as  considerably  simpler  than  solving  the 
RTE  in  its  original  form.  It  is  then  a  matter  of  recovering  the  radiance  field  via  an 
inverse  Fourier  transform  process: 

OO  OO  ^ 

/(f,  Co,  z)  =  dffiiK,  f],  z)  exp {i27r (k  -f+if-bj)}.  (42) 


2.6  Gaussian  Beam  Model 

The  method  proposed  for  solving  the  navigation  through  fog  problem  uses  a  laser 
source  operating  in  the  eye-safe  near-IR  spectral  region,  with  a  wavelength  A  of 
approximately  1.5  ym.  This  laser  source  will  emit  a  beam  focussed  approximately 
1  m  from  the  laser  exit  aperture  to  enable  a  wide  beam  divergence.  This  condition 
greatly  simplifies  the  mathematics  of  the  propagating  laser  beam  as  shown  in  the 
following  section. 


2.6.1  Complex  Beam  Wave 

Siegman  (1986)  provides  a  standard  numerical  model  for  describing  the  propa¬ 
gation  of  a  TEMoo  mode  Gaussian  beam.  This  form  should  be  equivalent  to  the 
vacuum  solution  of  the  RTE.  However,  due  to  differences  in  terminology  it  will  be 
necessary  to  translate  out  of  the  common  description  used  in  laser  literature  to  a 
form  that  is  compatible  with  the  radiance  description  commonly  used  in  radiative 
transfer. 

We  begin  by  describing  the  commonly  used  complex  beamwave  form.  Siegman 
describes  his  beamwave  in  terms  of  u,  the  beamwave  complex  amplitude: 


u{x,y,z) 


go 

\l  -K  Wo  q{z) 


exp 


-jkz-  jk 


(43) 


where  q{z)  is  the  complex  radius  of  curvature  for  the  Gaussian  beamwave,  wq 
is  the  beam  width  at  zero  range,  qo  is  the  complex  radius  at  zero  range,  k  is  the 
radiation  wavenumber  (27r/A),  and  j  =  r  =  \/x^"+^  as  before.  This 
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beamwave  form  assumes  that  the  minimum  beam  width  occurs  at  a  range  of  zero 
and  the  beamwave  travels  with  its  center  along  the  positive  2  axis. 

In  Siegman’s  derivation  the  paraxial  approximation  is  used.  While  recent  observa¬ 
tions  by  Zeng  et  al.  (1999)  indicate  some  correction  to  the  paraxial  approximation 
is  necessary,  we  can  avoid  these  corrections  by  the  choice  of  vertical  polarization 
of  the  emitted  beam  and  the  consideration  that  typical  off-axis  positions  of  interest 
will  lie  primarily  along  the  horizontal  axis. 

Of  considerable  interest  is  the  form  taken  by  q{z),  given  by  Siegman  as  a  function 
of  w{z),  the  beam  width  at  range  z,  R{z),  the  real-valued  radius  of  curvature  at  2, 
and  wavelength  A  through  the  expression, 


1  _  _ ■  A 

q{z)  R{z)  ^ 'nvP{zy 

where  the  radius  R  and  width  w  parameters  are  given  by, 

R{z)  =  z  +  —,  w^{z)  =  Wq 

z 


(44) 


(45) 


These  expressions,  in  turn,  hinge  on  definitions  of  zr,  termed  the  Rayleigh  range, 
and  Wo,  the  minimum  beam  width  which  is  assumed  to  occur  at  range  z  =  Q.  zr 
and  Wo  are  related  through  the  equation 


Zr  = 


A  ■ 


(46) 


Through  these  definitions,  the  equation  for  q{z)  can  be  written  as 

q{z)  =  z  +  qo  =  z  +  j  Zr.  (47) 


2.6.2  Radiance  Form  for  a  Diverging  Gaussian  Beam 


To  use  functions  u  and  q{z)  in  a  radiative  transfer  application  we  will  need  to  de¬ 
construct  these  functions  to  determine  their  relationship  to  the  radiance  function 
/  used  by  the  radiative  transfer  method.  To  begin,  we  multiply  the  complex  am¬ 
plitude  by  its  complex  conjugate.  This  produces  a  real-valued  flux  density  for  the 
propagating  beamwave: 


J(r,  z)  = 


2P 

TTW^ 


exp 


(48) 


where  we  use  w  =  w{z)  from  equation  (45).  Integration  of  X(r,  z)  across  the 
(x,  y)  plane  reveals  that  P  (power)  represents  the  energy  flowing  past  a  point  z 
through  the  (x,  y)  plane  per  unit  time  (given  in  Watts).  J(r,  z)  thus  has  units  of 
irradiance  (W/m^).  A  normalized  form  of  the  I  function  is  shown  in  figure  6, 
which  is  essentially  equivalent  to  Siegman’s  figure  17.2.  This  figure  includes  the 
dashed  outline  of  a  “top  hat”  beam.  The  top  hat  beam  has  the  property  of  constant 
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Figure  6.  Gaussian  beam  transverse  irradiance  relationship. 


irradiance  across  a  circular  beam  profile  of  radius  w/ \/2,  with  irradiance  equal  to 
the  maximum  of  the  related  Gaussian  beam,  and  exhibits  the  same  total  flux  as 
the  Gaussian  beam.  Hence,  the  top  hat  beam  is  sometimes  used  to  approximate 
properties  of  the  Gaussian  beam  yet  is  simpler  to  work  with  mathematically. 

The  J(r,  z)  form,  however,  only  provides  magnitude  information.  It  does  not  pro¬ 
vide  the  directional  information  needed  to  specify  u  functional  dependence.  To 
assess  this  additional  dependence,  the  complex  wavefront  curvature,  q{z),  equa¬ 
tion  must  be  analyzed.  A  key  part  of  this  analysis  involved  examination  of  Sieg- 
man’s  figure  17.1  (see  figure  7).  This  figure  appears  to  show  that  a  single  beam 
curvature,  R{z),  applies  for  every  point  in  each  plane  perpendicular  to  the  z  axis, 
where  R{z)  is  given  in  equation  (45). 

More  rigorously,  introducing  the  definition  of  q~^{z)  from  equation  (44)  into 
equation  (43)  results  in 


-/  \  /2  qo 

u(x,y,z)  —  \ - rTT 

^  ^  Uwoq{z) 


exp 


-jkz  -  jk 


2R{z)  w‘^{z)\ 


(49) 


This  form  reveals  a  separation  between  real  and  imaginary  structures  as  a  function 
of  off-axis  distance,  r.  The  exp(-r^/u;^)  term  accounts  for  the  off-axis  amplitude 
reduction,  while  exp{-jkr‘^/2R)  accounts  for  the  off-axis  phase  modulation.  This 
form  is  based  on  the  Fresnel  approximation,  utilizing  a  wavefront  curvature  of  R 
and  the  assumption  that  r  R. 

Given  this  constraint,  it  is  nevertheless  clear  that  use  of  the  phaseterm 
exp{-jkr'^/2R)  implies  that  at  any  given  point  (x,  y,  z)  the  wavefront  of  the  prop- 


17 


Figure  7.  Gaussian  beam  transverse  spread  {W)  and  curvature  (R)  as  functions  of  range  (Z). 


agating  beam  will  have  some  normal  vector  which  we  can  model  using  the  triplet 
{x/R{z),  y/R{z),  1}.  From  the  standpoint  of  adding  a  directional  dependence  to 
the  radiance,  this  normal  vector  leads  to  the  addition  of  a  two-dimensional  Dirac 
delta  function:  5  {u-  [r/R{z)]}  to  describe  the  directional  properties: 


^  r 

X 

y 

5 

r  ^(^)J 

=  5 

1 

5 

R{z) 

The  vacuum  form  for  the  radiance  of  a  propagating  Gaussian  beamwave 
should  then  be  written  as, 


0,  Z)  = 


2F 


TTW‘^{Z) 


exp 


-2r^ 


w^{Z) 


u  — 


R(Z) 


(51) 


where  we  again  mention  the  caveat  that  the  beam  is  rapidly  diverging  and  travel¬ 
ling  in  the  positive  2;  direction. 


2.7  Radiative  Transfer  Solution 

Equation  (40)  in  section  2.5  provided  a  Fourier  version  of  the  small  angle  approx¬ 
imation  of  the  RTE.  In  this  section  this  equation  is  solved  using  a  method  similar 
to  that  outlined  by  Smirnov  (1964).  Further,  the  solution  of  this  process  must  re¬ 
duce  to  the  equation  (51)  vacuum  Gaussian  beamwave  solution  in  the  absence  of 
aerosols  (cr  =  0); 

Equation  (40)  consists  of  a  series  of  linear  operators  applied  to  the  transformed 
radiance  function.  The  results  of  this  analysis  are  similar  to  those  of  Zardecki 
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and  Tam  (1982),  but  utilize  a  different  Fourier  transform  kernel,  modified  source 
radiation,  and  a  more  accurate  phase  function  description.  Also  here,  a  detailed 
description  of  the  mathematics  is  provided. 


2.7.1  Implicit  Solution  Method 


Smirnov’s  approach  involves  the  solution  of  a  single  linear  partial  differential 
equation  through  its  relationship  to  the  solution  of  a  set  of  linear  differential  equa¬ 
tions.  This  approach  begins  by  considering  a  set  of  variables:  xq,  Xi,  ...,Xn.  These 
variables  are  related  to  one  another  within  a  system  of  differential  equations: 


dxo  _  dxi  _  dx-i 

'X^~Tx~'X^ 


dXfi 


(52) 


The  Xi  are  functions  of  the  Xj  (Xi  =  Xi  {xq,  xi,  ...,  x„),  X2  = 
X2  (xo,  Xi,  ...,  x„),  etc.).  There  are  n  4-  1  equal  ratios  within  this  system  of 
relations. 

Within  this  system  there  are  effectively  n  dependent  variables  and  one  indepen¬ 
dent  variable,  though  the  form  used  for  the  relations  in  the  system  is  the  same  for 
all  n  -f-  1  variables.  Nonetheless,  one  of  these  variables  could  be  chosen  as  the 
independent  variable.  For  example,  if  the  rco  variable  is  chosen  to  be  the  indepen¬ 
dent  variable,  then  the  other  variables  can  be  written  as  xi(xo),  X2(xo), ...,  x„(xo). 
If  we  are  able  to  write  out  explicit  equations  for  each  of  the  dependent  variables 
as  functions  of  the  single  independent  variable  then  the  system  of  equations  (52) 
is  solved. 

Smirnov  posits  such  a  system  of  equations  and  its  solution  and  then  makes  a  sec¬ 
ond  point:  He  assumes  that  given  such  a  system  of  equations  one  can  produce 
a  function,  call  it  E{xo,  xi,  ...,  Xn).  This  function  has  the  property  that  it  is  not 
identically  zero.  Yet,  if  any  solution  to  the  system  of  equations  (52)  is  substituted 
for  all  of  the  dependent  variables  in  that  system,  then  E  evaluates  to  a  constant. 
Smirnov  calls  such  a  function  an  integral  of  the  system  of  equations. 

Viewed  another  way,  the  function  E  can  be  considered  to  define  a  hypersurface 
within  the  hyperspace  (rro,  xi,  X2,  ....  Xn).  This  hypersurface  corresponds  to  a 
solution  of  the  system  of  equations  (52)  that  satisfies  the  condition. 


^(xqj  Xi,  X2,  Xyj)  C* 


(53) 


In  the  analysis  to  follow,  Smirnov  exploits  the  properties  of  this  notional  function, 
E.  However,  given  the  nature  of  this  analysis,  he  never  needs  to  actually  solve 
explicitly  for  function  E  in  terms  of  the  set  of  variables  xo, ...,  Xn.  The  first  step  in 
this  exploitation  process  involves  generating  the  total  derivative  of  c,  with  respect 
to  the  chosen  independent  variable  xq. 


dE  dE  dE  dxi  dE  d^  dE  dxn 

dxQ~  dxo  dxidxo  8x2  dxo  *  dxo 


(54) 
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Or  in  condensed  notation, 


cE  _  A  52 
dxo  ^  dxi  dxo 


(55) 


In  this  equation  the  terms  dxi/dxo  merely  specify  the  rate  of  change  in  one  of  the 
variables  (xi)  relative  to  the  change  in  the  variable  xq.  One  may  imagine  some 
random  positional  displacement,  within  the  hyperspace.  Essentially  this  dis¬ 
placement  is  completely  arbitrary,  and  hence  the  change  in  E  may  also  be  nonzero. 
However,  if  the  direction  of  displacement  is  determined  by  calculating  the  quan¬ 
tities  dXi/dxo  using  solutions  of  the  system  of  equations  for  the  Xj  in  terms  of  xo, 
then  the  net  change  in  5  must  equal  zero!  Multiplying  both  sides  of  equation  (55) 
by  dxo  we  generate, 

where  the  use  of  the  subscript  S  implies  that  the  dependent  dxi  are  constrained  to 
vary  only  along  the  solution  hypersurface. 

To  exploit  the  development  in  equation  (56)  an  important  auxiliary  development  is 
now  required.  In  this  development  we  assign  a  quantity,  dX,  to  all  the  equal  ratios 
given  in  equation  (52): 

dxQ  dx\  dxi  dXfi 

dX  -yr  y  y  •’*  V*  ^  ^ 

Ao  Ai  Ai 

Since  dX  is  a  nonzero  quantity,  it  can  be  divided  into  equation  (56),  and  terms  can 
be  cancelled,  to  obtain. 


(56) 


i=0  1=0 


(58) 


We  now  compare  the  appearance  of  equation  (58)  to  equation  (40),  the  Fourier 
version  of  the  RTE.  The  form  of  equation  (58)  is  nearly  identical  to  the  equation 
to  be  solved  except  that  the  partial  derivatives  in  equation  (58)  only  appear  as 
operators  on  the  integral  function  S,  whereas  the  equation  (40)  has  the  form. 


5x„  5Xn  dXn 


+  Fn  =  0, 


(59) 


where  the  Yi  are  functions  of  xo,  xi  ,  ...  ,  and  x„,  and  the  partial  differentials 
operate  on  the  variable  Xn,  which  itself  appears  in  the  Yi  functions.  By  contrast,  the 
Xi  functions  do  not  depend  on  E.  It  is  therefore  necessary  to  relate  equation  (59) 
to  the  form  of  equation  (58)  by  replacing  the  partial  derivatives  of  x„  with  partials 
involving  E. 


To  perform  this  transformation,  equation  (56)  is  again  examined  under  the  condi¬ 
tion  of  holding  all  but  two  variables  constant.  Let  the  two  free  variables  be  chosen 
as  Xi  and  x„.  Therefore,  from  equation  (56), 


dE  ,  9E 

- —  dxn  +  —dxi  =  0. 

UXfi  OXi 


(60) 
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or, 


(61) 


dxn  _  dEfdxi 
dxi  dE/dxn 

In  this  case  only  Xn  and  Xi  are  allowed  to  vary  and  thus  this  result  is  a  partial 
derivative  (dxn/dxi).  Smirnov  calls  this  the  rule  for  differentiation  of  implicit 
functions.  This  rule  allows  the  partial  derivatives  of  Xn  to  be  replaced  by  partials 
with  respect  to  E.  Equation  (59)  is  thus  transformed  into. 


°  dxo  ^  dxi 


o«— > 

+  y„^  =  o. 

OXn 


(62) 


We  thus  have  a  method  of  relating  equation  (59),  which  has  the  same  form  as  the 
differential  equation  we  wish  to  solve,  to  equation  (58),  an  equation  which  we 
have  causally  connected  to  the  system  of  equations  (52).  In  particular,  we  have 
established  connections  between  — y  Xi,  for  all  i  <  n\  and  Yn  Xn-  Any 
partial  differential  equation  that  can  be  written  as  a  system  of  linear  operators  in 
the  form  of  equation  (59)  may  be  related  to  a  system  of  equations  in  the  form  of 
equation  (52).  The  expectation  is  then  that  this  related  system  of  equations  will  be 
simpler  to  solve  than  the  original  partial  differential  equation. 

When  we  apply  this  series  of  steps  to  our  problem  the  first  step  is  to  assign  2:  as 
the  independent  variable  (xq).  The  remaining  variables  are  assigned  as,  xi  =  r]^, 

X2  =  riy,  and  xz  =  Xn  =  I.  Next,  Yi  functions  are  assigned: 


-Xo  =  yo  =  1; 

Xl  =  Y\  =  ^35) 


X2  — 

X3  =  ¥3=  [a-  waPf\  i. 


(63) 


Introducing  expressions  for  the  Xfs  into  the  related  system  of  equations  generates 
three  independent  relations,  the  first  two  of  which  are 


-dz  =  ^ 


(64) 


Integration  of  these  first  two  functions  produces  a  relationship  between  the  angular 
frequency  vector,  ff,  the  spatial  frequency  vector,  k,  and  range  2:: 

f]{z)  =  C,j  -  Kz,  (65) 

where  (7,,  is  a  constant  of  the  integration. 

This  result  appears  to  indicate  that  77  is  a  completely  dependent  function  of  k  and 
z,  but  this  is  not  the  case.  Rather,  if  is  only  dependent  on  z  and  k  in  terms  of  its 
positional  variability.  In  fact,  Cj,  may  freely  vary  as  long  as  it  is  independent  of  2:. 
In  particular, 

rf(0)  =  Cr,  =  rfo.  (66) 
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Essentially  the  restrictions  on  ff  imply  that  propagation  under  the  radiative  transfer 
model  involves  a  corpuscular-type  assumption  —  that  once  a  photon  is  launched  in 
a  given  direction  it  does  not  arbitrarily  begin  propagating  in  some  other  direction 
as  it  travels  forward  in  z.  Thus,  diffraction  effects  are  not  treated  directly.  This  is 
a  limitation  of  the  radiative  transfer  method  which  ignores  large-scale  diffraction 
effects.  However,  small-scale  diffraction  is  modeled  indirectly  through  use  of  the 
scattering  phase  function. 

A  second  choice  for  defining  Cj,  is  at  path  position  z  =  Z.  At  this  range, 

ff{Z)  =ffz  =  ffo-KZ,  (67) 

Cj,  =  ffo  =  ffz  +  Z.  (68) 

These  two  means  of  expressing  Crj  are  significant  because  we  will  deal  with  spe¬ 
cific  examples  of  propagation  where  we  are  interested  in  the  propagated  radiance 
field  at  ranges  0  and  Z,  and  defining  Cr,  in  terms  of  these  ranges  is  important  in 
understanding  the  nature  of  the  propagating  radiance  relative  to  these  particular 
ranges.  Once  C,,  is  determined,  equation  (65)  indicates  ffs  range  dependence  at 
all  other  ranges.  In  the  forward  propagating  laser  beam  case,  ffz  is  used  as  the  an¬ 
gular  frequency  variable  since  the  properties  of  the  propagated  signal  as  it  enters 
the  retroreflector  will  be  of  most  interest. 

Having  thus  chosen  the  rfz  representation  for  Cr,,  fj{z)  is  written  as, 

ff{z)  =  f]z  +  K{Z-z).  (69) 


2.7.2  Fourier  Solutions 

The  third  and  final  relationship  derived  from  equation  set  (63)  is  given  by, 

^  = —a{l-wPf)dz,  (70) 

i 

which  has  the  solution, 

I  =  Iq  exp  1“  /  ^ 

Note  that,  similar  to  the  case  with  ff,  Iq  need  only  be  a  constant  with  respect  to  the 
z  variable.  This  is  illustrated  as  we  consider  several  solutions  to  this  equation  for 
various  propagation  conditions. 

2.7.2. 1  Vacuum  Solution 

As  observed  previously,  the  formula  given  for  in  equation  (51)  must  be  the 
solution  of  equation  (71)  in  the  case  of  a  vacuum  atmosphere  (i.e.,  when  a  =  0). 
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But  when  a  =  0  the  argument  of  the  exponential  in  equation  (71)  equals  zero.  In 
this  case  I  =  /q  in  the  Fourier  representation  of  the  solution.  Thus,  the  Fourier 

transform  of  must  be  Iq,  a  constant  with  respect  to  ^  according  to  the  radiative 
transfer  method. 

Rewriting  the  vacuum  solution  at  range  z  =  Z, 


I^if,  (I!,  Z)  = 


2P 


-Kw\Z) 

this  function  Fourier  transforms  to, 


exp 


■  -27-2  ■ 

X 

w^{Z) 

0 

[  m\ 

(72) 


fi,Z)^P  exp 


^2  ^2  ^2 
2 


(73) 


where  A  is  the  magnitude  of  the  vector, 

K  =  K^-fi/R.  (74) 


Yet  both  R  (the  radius  of  wavefront  curvature  from  equation  (45))  and  fj  are  func¬ 
tions  of  position,  2:,  along  the  path.  At  first  glance,  then,  appears  to  violate 
the  requirement  that  its  transform  be  constant  with  respect  to  range.  This  diffi¬ 
culty  arises  because  the  actual  equation  for  the  propagating  beamwave  does  not 
truly  follow  radiative  transfer  theory’s  corpuscular  propagation  assumption.  While 
this  difference  cannot  be  completely  resolved,  wA  is  approximately  constant  with 
range  as  long  as  zr  z.  Considered  another  way,  as  long  as  ttuiq  <C  Xz,  such 
that  the  characteristic  width  of  the  beam  at  its  focal  point  is  small  compared  to  a 
Fresnel  zone  length,  then  the  beam  can  be  considered  to  be  nearly  a  point  source. 
This  condition  is  possible  for  the  laser  system  contemplated  since  the  beam  is 
focused  at  a  short  range  from  the  laser  aperture. 


Under  these  conditions. 


wlz)  =  Wo  \  1  -1-  -5- 

V  4 


Wo 


where  wo  varies  as. 


Wo  = 


Xzr 


Wq 

zr 


zr 


A 


TT  Zr  TT  Wo 

Conversely,  A  varies  approximately  inversely  with  range, 
+  k{Z-z)] 


A  —  K.  -\- 


^  ,  Vz  ,  _  {fjz  +  Z k)  ffo 

K-\ - K - K  =  -  =  — . 

Z  2:  2:  2: 


(75) 


(76) 


(77) 


^  [1  +  (4A^)] 

Combining  the  influences  of  w  and  A,  the  argument  to  the  exponential  in  the 
transformed  vacuum  radiance  is  given  by, 


TT 


—  W^  A 


,2  A  2 


(78) 


2  “  ** 

which  is  independent  of  range  2:.  Thus,  under  the  limitations  mentioned,  the  vac¬ 
uum  solution  of  equation  (72)  is  suitable  as  a  solution  to  the  RTF  in  Fourier  space. 
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2.112  Unscattered  Solution 


Once  the  vacuum  form  of  the  solution  is  known,  the  unscattered  solution  may  be 
immediately  written  as, 

/^(f,  dJ,  Z)  =  I^{f,  iv,  Z)  exp  [-r(Z)] ,  (79) 

where 

z 

T{Z)  =  Jadz,  (80) 

0 

is  the  optical  depth  encountered  by  a  beam  travelling  through  distance  Z. 

In  Fourier  space  this  equation  is  written  as, 

F  =  exp[-T{Z)],  F  =  io.  (81) 

This  result  is  obtained  by  setting  the  single  scattering  albedo,  zu,  to  zero. 


2.73  Scattered  Solution 


The  solution  for  the  scattered  radiation  case  is  considerably  more  complex.  We 
begin  by  resolving  it  into  a  summation  of  scattered  and  unscattered  components: 


I  = 


F  +  F] 

tu  ,  w 


=  F+F  [5-1]; 


=  F  +  F 


exp  ^  J  waPf  dz 
0 


-  1 


(82) 


where  F  is  the  scattered  solution  and  S  characterizes  the  spatial  and  angular 
frequency  effects  of  the  scattering. 

It  will  be  necessary  to  evaluate  S  in  order  to  determine  its  influence  on  the  radi¬ 
ance  field  arriving  at  the  retroreflector.  To  simplify  this  task,  thereby  reducing  the 
number  of  free  parameters  involved  in  this  study,  it  is  assumed  that  the  scenario 
studies  will  be  one  involving  an  aircraft  landing  at  an  airfield  enshrouded  in  fog. 
As  the  aircraft  approaches  the  field,  at  some  point  the  plane  will  be  above  the  fog. 
Then,  as  the  plane  makes  its  final  approach,  it  becomes  completely  immersed  in 
the  fog.  Thus  the  optical  path  may  be  modeled  to  first  order  as  consisting  of  an 
initial  portion  free  of  fog  and  a  latter  portion  (containing  the  retroreflector)  im¬ 
mersed  in  fog.  The  fraction  of  the  modeled  volume  that  is  free  of  fog  will  range 
from  zero  (completely  free  of  fog)  to  unity  (completely  immersed  in  fog).  For 
simplicity,  the  volume  is  assumed  to  have  uniform  density  and  scattering  proper¬ 
ties  throughout  the  fog  region.  Zf  is  designated  as  the  distance  from  the  laser  focal 
point  to  the  forward  edge  of  the  fog  enshrouded  portion  of  the  modeled  volume. 
The  dimensionless  quantity  Q  =  {Z  -  Zf)/Z  measures  the  fraction  of  the  optical 
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path  containing  fog.  The  fog  aerosol  properties  are  assumed  constant  in  both  x 
and  y,  as  characterized  by  parameters  a,  w,  and  the  coefficients  constituting  the 
phase  function  description.  The  remainder  of  the  optical  path  is  assumed  to  be 
filled  with  an  optically  thin  aerosol  whose  properties  represent  a  virtual  vacuum. 

Given  these  limiting  approximations,  the  scattering  solution  reduces  to  an  evalua¬ 
tion  of  the  function  S  and  the  inverse  transform  of  the  resulting  function.  Because 
the  forward-scattering  phase  function  was  assumed  to  be  a  sum  of  Gaussian  terms, 
it  is  possible  to  exactly  represent  S. 


2.1.4  Error  Function  Application 


In  solving  for  <S,  integrations  involving  Gaussian  functions  must  be  performed. 
These  integrals  use  the  error  function,  given  by 


2  ^ 

erf(x)  =  -j=  j  exp  (— dt. 

Defining  a  Gaussian  function  according  to  Gaskill’s  (1978)  definition, 

Gaus(x)  =  exp  (-tt  , 
the  integral  of  the  Gaussian  can  be  written  as, 

T  (X2,  Xi)  =  j  Gaus(w)  du  =  -  [erf  (V7rx2)  -  erf  (\/^Xi)] . 

Xi 

Integrating  along  the  path  for  the  scattering  effects  in  S, 


(83) 


(84) 


S  =  exp  < 


Z 

wa  j  Pf  {\ffz  +  k[Z-  z]\)  dz 


M 


=  n^xp^ 


Zf 

i  Am  ■CUT 


m=l  I  QQv 


exp  [-Trp^  (l  -  T  [Qq^  -f  ppjX  ,  (85) 


where  r  =  aQ  Z,qm  =  ZKfam,  Pm  =  Vz/otm,  //=(«•  rfz)  /  {ktjz)  is  the  cosine 
of  the  angle  between  the  k  and  ffz  vectors,  and  k  and  rjz  are  magnitudes  of  their 
respective  vectors. 

Using  this  result,  the  solution  in  Fourier  space  may  be  written  as 


I  =  P  exp  [-Ag  (9^  -I-  2pqp  -f  exp(-r)  x 


n  ^  [QQm  +  PPm,  PPm]  exp  [“TTp^  (l  -  |  , 


m=l 


where 


Ag  = 


2wl' 


(86) 


(87) 


P,  Aq,  Q,  t,  and  zu  are  constants  of  the  laser  propagation  scenario,  while  q,  p,  and 
p  provide  the  variability  associated  with  k,  and  ffz  in  terms  of  unitless  parameters. 
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2.7.5  The  Inverse  Transform 

Because  the  Fourier  representation  of  the  radiance  solution  has  factors  whose 
magnitudes  peak  at  different  points  in  frequency  space,  numerical  methods  will 
be  necessary  to  evaluate  the  inverse-transformed  radiance  field.  In  its  full  form, 
the  inverse  transform  is  given  by 

00  OO  ^ 

J(f,  u,Z)  =  /(k,  nz^  Z)  exp  0'27r  (k  -f+m-  •  (88) 

— OO  —CX) 

However,  processing  this  transform  can  be  simplified  by  reducing  the  order  of 
the  integral  by  one.  This  step  is  possible  because  equation  (86)  only  involves  the 
quantities  p,  q,  and  /x.  Here,  let  /x  =  cos(<p)  define  the  angle,  <p  formed  between 
the  vectors  k  and  ff.  Given  the  symmetries  of  the  problem,  it  is  possible  to  solve 

for  the  radiance  field  in  terms  of  r  =  |r1,u;  =  Iwl,  and  the  angle  formed  between  r 

and  w  (call  it  a),  instead  of  f  and  xD  in  fully  four-dimensional  form.  This  reduction 
acknowledges  the  radial  symmetries  that  exist  at  off-axis  positions. 

To  integrate  out  one  of  these  four  dimensions,  consider  figure  8.  From  the  ge¬ 
ometry  in  this  figure,  the  dot  products  contained  in  the  complex  argument  of  the 
exponential  in  equation  (88)  can  be  written, 

K  .  r  =  K r  cos(a  +  c  —  (p  —  ip),  (89) 

jfz  •  w  =  pz  w  cos(c  -  Ip) .  (90) 


The  form  that  these  dot  products  take  suggests  that  if  the  multiple  integral  in 
equation  (88)  were  expressed  in  polar  form  that  the  ip  variability  in  the  integral 
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can  be  immediately  integrated  out.  In  polar  form  the  inverse  transform  becomes 


oo  r  oo  r  27r  r  27r 


//<iic//di;zi...l-/ 


d(p 


7]  dr) 


ndn,  (91) 


-oo  -oo 


this  fixed  structure  ttte  start  and  finish  points  of  the 

TO  accomplish  the  integmtion.  the  flmt 

We^X  ifaltonirtfie 

P  is  introduced,  where  ^  c  _  ip  - 1/^ 

Since  ^  is  constant  while  ^  is  varild.  0  =  -d^-  Tlterefom,  upon  replacing  ^ 
with  p,  the  innermost  integral  becomes 

J exp  {j  27r  [k  r  cos{p)  +  r?2  w  cos{p  +  (^  -  «)]} 

This  integral  I  solved  analytically  using  a  result  from  Goodman  (1968): 

j  exp  {-j  L  cos{9'  -  <!>')]  d9'  =  2tx  Jo{L)  ,  (94) 

where  Jo  is  a  zero  order  Bessel  function  of  the  first  kind.  Integral  (93)  can  be 
related  to  integral  (94)  by  observing  that  the  quantity 

_  —  L  cos(0'  -  (j>') 

27r 

is  simply  the  ^-component  of  a  vector, 


—  (-L  cos  [-</)'] ,  -L  [-sin0']) 
27r 


(95) 


that  is  being  rotated  about  the  origin,  with  an  offset  angle  P'. 

S(£odmrs'iuor?hrge™^^^^ 

?he  vector  U  in  that  figum  is  composed  of  the  vector  sum  of 


=  (/cr,  0) , 


(96) 
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Figure  9.  Orientation  of  vector  U  relative  to  component  vectors  KV  and 

and 

Tjz(j^  =  cos(<^  —  a),  sin(^  —  a)) .  (97) 

Goodman’s  integral  result  is  independent  of  the  value  of  4>'  since  a  unit  circle  is 
being  traversed.  Similarly,  U  ignores  any  initial  offset  between  the  x  axis  and  the 
/cr  vector.  The  only  quantity  of  interest  is  the  magnitude  of  U,  given  by 

U{Kr,  Tjztj,  ^  —  a)  =  1^1  =  +  2 nr Tjzu)  cos{(p  —  a)  +7/f  (98) 

Hence,  the  innermost  integration  of  equation  (91)  evaluates  to 

27r 

J exp  {j  27r  [fcr  003(7)  +'nz‘^  003(7  +  <P  —  o)]}  <^7  =  27r  Jo{2ir  U).  (99) 

0 

Eliminating  this  inner  integral  significantly  simplifies  the  calculations  necessary 
to  evaluate  the  radiance  field  for  a  specific  (r,  u,  a)  triplet.  Simultaneously,  in¬ 
tegration  with  Ip  eliminates  the  need  to  work  with  a  complex  exponential  kernel 
in  the  inverse  transform.  The  resulting  Bessel  function  ensures  that  the  resulting 
integration  will  be  purely  real-valued. 

This  result  also  has  immediate  implications  for  the  resulting  inverse  transform 
when  both  r  and  w  are  zero.  For  this  case  the  Jo  term  evaluates  to  unity  and  all 
dependence  on  the  angle  a  drops  out  of  the  resulting  inverse  transform. 

2.8  Numerical  Integration 

Having  removed  one  order  of  integration,  we  must  now  evaluate  the  triple  integral 
involving  k,  rjz,  and  <p. 
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To  accomplish  this  task  various  numerical  methods  were  considered,  such  as  those 
contained  in  Press  et  al.  (1992).  However,  of  the  available  sources  of  routines  none 
were  found  to  be  applicable.  Hence,  this  section  describes  a  new  technique  for 
treating  this  task. 

First,  as  previously  stated,  different  terms  in  the  integration  kernel  have  peaks 
at  different  locations  in  the  integration  volume.  For  instance,  the  Jo  term  peaks 
where  C/(/cr,  ^  -  o)  =  0.  This  condition  occurs  along  the  plane  kt  =  rizoJ 

since  there  will  always  be  a  v?  value  for  which  cos{(p  —  a)  =  —  1.  Hence,  following 
integration  by  (p,  there  will  be  a  resonance  along  nr  =  rjzu}-  (This  resonance  also 
encompasses  the  specific  condition:  nr  =  rjzcj  =  0.) 

On  the  other  hand,  the  term 


exp 


2wo 


{Zn  +  rizf 


(100) 


has  a  resonance  when  ZK  =  r)z. 

Comparing  these  two  expressions,  there  is  a  double  resonance  for 


Not  surprisingly,  this  peak  describes  the  direction  of  geometrical  propagation  for 
imscattered  radiation. 

The  remaining  S  term  in  equation  (86)  provides  the  scattering  dependence  for  I. 
This  function  has  a  very  complex  behavior,  depending  on  the  particular  combi¬ 
nation  of  parameters  describing  the  aerosol  and  path.  Basically,  all  the  features 
describing  the  optical  path  are  here:  the  thickness  of  the  aerosol,  the  scattering 
properties,  the  portion  of  the  path  in  the  fog,  and  the  angular  and  positional  struc¬ 
ture. 

Even  so,  a  few  generalizations  are  possible.  For  example,  when  /x  =  —  1,  the  T 
function  tends  to  exhibit  a  broad  peak.  Also,  when  k  — ^  0, 


QQm 


exp[-7r/x2p^  . 


(102) 


Thus,  some  cancellation  occurs  with  the  exp  [— (1  -  /x^)]  term,  resulting  in  a 
exp  (— TTjo^)  =  exp  (— TTT/l/a^)  dependence  around  q  =  0.  Unfortunately,  these 
general  observations  do  not  assist  in  the  particular  evaluation  of  the  Fourier  ker¬ 
nel.  Rather,  a  triple  integral  must  be  evaluated  numerically,  including  k  and  rjz 
integrations  out  to  infinity. 

Press  et  al.  (1992)  discuss  these  types  of  improper  integrals  in  their  numerical 
recipes  book.  Their  method  proposes  dividing  the  integral  into  two  parts,  delimited 
by  a  midpoint.  In  the  lower  section  the  integral  is  evaluated  normally.  In  the  upper 
section  a  change  of  variables  is  used,  wherein  the  inverse  of  the  original  variable  is 
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introduced.  Through  this  change,  the  limits  of  integration  are  modified  to  a  finite 
interval.  Mathematically,  this  is  written 

OO  Xm  OO 

j  f{x)dx  =  j  f{x)dx+  j  f{x)dx 

0  0  ^rn 

Xm 

f{x)dx+  j  f{y~^)y~‘^dy  (103) 

0  0 

where  y  =  x~^.  But  this  analysis  only  considered  one-dimensional  integration. 

Since  the  triple  integral  of  J  involves  integration  to  infinity  in  two  dimensions, 
to  use  separate  divisions  of  rjz  and  k  space  to  perform  this  integration  would 
result  in  four  different  integration  regions.  The  problem  with  this  approach  arises 
in  terms  of  the  resonance  regions  in  the  integrand.  Because  these  regions  tend 
to  emanate  from  the  origin,  they  become  difficult  to  resolve  accurately  at  high 
argument  values.  At  distances  further  from  the  origin  the  nonzero  region  tends 
to  become  closer  to  the  resonant  line  itself,  meaning  that  only  a  small  portion  of 
the  overall  integration  region  is  actually  responsible  for  the  integral  results.  The 
remainder  of  the  area  has  an  effective  integrand  value  of  zero.  Using  rectangular 
integration  regions,  much  larger  numbers  of  samples  are  required  to  accurately 
capture  the  behavior  of  the  function  at  high  argument  values. 

To  overcome  this  difficulty,  variables  p  and  q,  which  already  represent  polar  mag¬ 
nitudes  of  T)z  and  k  variables,  are  combined  into  yet  another  polar  variable,  p.  To 
permit  this  representation  it  first  seemed  logical  to  modify  the  integration  space 
from  the  rjz  and  k  variables  to  the  dimensionless  p  and  q  variables: 

OO  OO 

jj fi  («,  Vz)  « dK T]z  drjz  =  Z-^  Jj fi  p'j  pdpq  dq.  (104) 

-OO  -OO 

Introducing  Wr  =  r/Z^  the  arguments  of  U  become 

U(Kr,  rjzijJi  <i>  —  a)  =  U{qu)r,  pu),  (j)  —  a).  (105) 


Thus, 


OO  OO  zn  ^ 

=  ^  J QdqJ pdp  J d(t>I\p,  q,  cos(0)]  Jo  [27rl7] .  (106) 


Introducing  polar  coordinates  p  and  6,  such  that  p  =  p  cos{9),  q  =  p  sin(0),  we 
can  rewrite  the  integration  with  respect  to  p  and  q  as. 


OO  OO  OO  7r/2 

Jqdqjpdp.,.  =  J P^ dp  J  - sin(20) d0...  (107) 

0  0  0  0 

Having  obtained  this  result,  we  can  now  make  a  few  observations.  First,  let  us  con¬ 
sider  the  integration  kernel  on  the  RHS  of  equation  (105).  Due  to  the  p^  weighting 
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we  do  not  expect  to  find  the  most  significant  wavenumber  response  near  the  origin 
nor  along  either  the  +p  or  +q  axes.  Rather,  we  expect  the  most  significant  region 
to  lie  along  the  diagonal  p  =  q-  This  consideration  will  to  some  degree  dictate  the 
means  used  as  we  attempt  to  numerically  evaluate  the  resulting  integral.  Second, 
the  use  of  the  q  variable  implies  a  similarity  of  results  measured  at  different  scales. 
Since  off-axis  variables  are  all  expressed  relative  to  the  overall  path  length,  we  can 
study  the  effects  of  scattering  of  fogs  in  the  laboratory  and  rescale  the  lengths  to 
apply  the  measured  results  to  large  scale  problems. 

The  computer  procedure  for  evaluating  the  remaining  (p,  6)  representation  of  the 
integral  divides  the  quarter  plane  integration  region  into  sections  according  to  fig¬ 
ure  10.  The  first  region,  labeled  “1”  in  this  figure,  covers  the  full  range  of  6  vari¬ 
ation.  The  limit  of  p  for  this  region  is  chosen  such  that  I  has  been  significantly 
reduced  below  its  value  at  p  =  0.  However,  because  of  the  p®  dependence,  the 
remainder  of  the  quadrant  must  also  be  evaluated  using  type  2  and  type  3  regions 
as  labeled  in  figure  10.  The  latter  region  types  are  created  by  subdividing  the  6 
integral  into  a  series  of  segments,  A6.  In  each  of  these  segments  successive  in¬ 
tegration  regions  of  type  2  are  evaluated,  beginning  at  the  boundary  of  the  type 
1  integration  region,  and  continuing  outward  until  the  segment  integration  result 
falls  below  some  threshold.  Once  this  threshold  is  crossed,  the  remainder  of  the 
Ad  interval  can  be  evaluated,  fi-om  pm  (the  maximum  p  value  for  this  A6  seg¬ 
ment)  out  to  infinity,  via  introduction  of  the  variable  ^  =  p~^.  This  final  integra¬ 
tion  avoids  the  problems  associated  with  integrating  over  both  rjz  and  k  because  a 
single  diagonal  region  is  considered. 


Figure  10.  Division  of  the  quarter  plane  in  J)  and  Q,  using  regions  delineated  by  constant  surfaces  in  the 
cylindrical  coordinates  p  and  9. 


The  practical  reason  for  dividing  the  integration  task  between  a  single  type  1  re¬ 
gion  and  many  type  2  regions  follows  from  the  nature  of  the  integrand  being  eval- 
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uated.  Because  many  of  the  diagonals  will  make  little  or  no  contribution  to  the 
overall  integral,  it  is  possible  along  many  of  these  A6  regions  to  quickly  move 
from  making  type  2  calculations  to  making  a  single  type  3  calculation.  However, 
there  will  probably  be  one  or  two  diagonals  that  have  significant  contributions 
out  to  large  p  values.  Then,  most  of  the  computation  time  will  be  focused  on  the 
few  diagonals  where  more  of  the  non-negligible  integrand  lies.  This  optimizes  the 
computational  accuracy  and  speed  of  the  technique. 


3.  Retroreflector  Interactions 


In  the  previous  chapter  the  propagation  of  forward  scattered  radiation  to  a  retrore¬ 
flector  was  considered.  In  this  chapter  the  interaction  between  the  incident  laser 
beam  and  the  retroreflector  is  considered.  In  the  analysis  of  this  interaction  a  spe¬ 
cific  retroreflective  material,  similar  to  the  materials  used  in  roadside  signs,  is  con¬ 
sidered.  This  material  is  composed  of  a  series  of  miniature  retroreflectors  (mini¬ 
reflectors)  embedded  in  a  plastic  sheathing.  This  retroreflective  material  can  be 
used  in  a  standalone  role  or  may  be  fashioned  into  a  larger  retroreflector  that  com¬ 
bines  the  microproperties  of  the  retroreflective  material  with  the  properties  of  a 
macro  retroreflector. 

At  the  material  level  the  retroreflective  material  (retro-material)  is  composed  of 
a  hexagonal  grid  structure  of  adjacent  mini-reflectors  similar  to  that  shown  in 
figure  11. 


Honeycombed 

Retroreflectors 


Interreflector 

Gaps 

Figure  11.  Set  of  gridded  mini-reflectors  as  they  appear  on  reflective  material. 

In  this  material  the  mini-reflectors  are  embedded  within  a  layer  of  plastic  that 
coats  the  entire  surface  of  the  material.  The  mini-reflectors  appear  to  be  circular 
in  cross-section  and  have  some  separation  distance  between  adjacent  retroreflect- 
ing  regions.  Thus,  there  is  some  fraction  of  the  reflective  material  surface  over 
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which  the  incident  signal  will  not  be  retroreflected.  Tests  of  retro  material  sam¬ 
ples  also  indicate  that  there  is  a  significant  specularly  reflected  component  when 
the  material  is  illuminated  with  a  handheld  laser  pointer.  Lastly,  there  is  a  sub¬ 
stantial  retroreflective  component,  but  it  is  apparent  that  the  combined  effects  of 
the  plastic  sheathing  and  the  inter-reflector  gaps  may  reduce  the  total  amoimt  of 
energy  retroreflected  from  the  material  to  below  50  percent. 

The  properties  of  the  retroreflective  regions  can  be  understood  by  assuming  that 
each  individual  mini-reflecting  region  is  a  circular  entrance  aperture  opening  onto 
three  orthogonally  oriented  mirrored  surfaces.  The  geometry  of  one  of  these  re¬ 
flectors  is  seen  in  figure  12.  This  figure  shows  views  of  the  three  orthogonal  mir¬ 
rored  surfaces  fi'om  both  the  front  and  the  side. 


Figure  12.  Cleometry  of  a  single  mini-reflector. 


In  the  figure,  variables  q',  r',  and  s'  are  introduced:  r'  is  the  radius  of  an  individ¬ 
ual  mini-reflector  as  measured  in  the  main  surface  plane  of  the  retro-material;  s' 
represents  the  depth  of  the  mini-reflector;  and  q'  measures  the  diagonal  distance 
between  points  P  and  A,  representing  the  seam  between  adjacent  mirrored  sections 
of  the  mini-reflector. 

Variables  g'  and  s'  are  related  to  r'  through  the  relations,  q'  =  r'  and  s'  = 

r'  J\/2.  These  relations  are  found  by  computing  the  distance  between  points  A 
and  B.  The  angle  AOB  has  an  arc  measuring  120°.  The  two  legs,  OA  and  OB, 
in  triangle  AOB  are  both  of  length  r'.  Thus  the  distance  AB  can  be  computed  as 
^  r'.  The  right  triangle  APB  has  two  legs  AP  and  BP  of  length  q'  and  the  third  leg 
of  length  AB.  The  q'  value  immediately  follows,  and  s'  can  be  evaluated  through 
computations  made  on  triangle  AOP. 

While  complicated  analyses  would  likely  produce  the  same  result,  simple  peering 
through  a  large  scale  retroreflector  reveals  that  beyond  the  reflection  point,  one 
sees  an  image  of  the  entrance  aperture  a  distance  2s'  “into”  the  retroreflector. 
Figure  13  illustrates  an  equivalent  model  of  this  situation.  The  fraction  of  light  that 
exits  the  system  following  passage  through  the  entrance  pupil  at  a  given  angle  of 
incidence  is  modeled  as  though  the  light  were  passing  through  a  circular  aperture 
of  radius  r',  through  a  short  tunnel  of  length  2s',  and  out  a  second  circular  aperture 


Figure  13.  Effective  geometry  of  a  retroreflector  considering  the  fraction  of  incident  energy  that  returns 
in  the  retro-direction. 


of  radius  r'.  In  effect,  the  mirror  properties  of  the  retroreflector  produce  an  image 
of  the  entrance  aperture  a  distance  s'  beyond  the  deepest  point  (P)  within  the 
retroreflector.  There  will  thus  also  be  a  maximum  angle  for  which  a  portion  of 
the  light  can  enter  the  system  at  the  edge  of  the  entrance  aperture  and  reach  the 
center  of  the  retroreflector.  All  light  entering  at  a  greater  angles  of  incidence  will 
not  strike  all  three  reflective  surfaces  and  cannot  retroreflect. 

Mathematically,  the  portion  of  the  energy  which  returns  in  the  retro-direction  is 
given  as  the  fractional  overlap  area  between  two  circles  of  equal  diameter.  This 
fraction  is  given  by  the  formula, 

fv{u)  —  —  [cos“^(m)  —uVl  —  cyl(«/2),  (108) 

Tj-  L  J 


where  fv  is  the  vignetting  factor,  u  =  tan(0')/x/2,  and  cyl(u)  is  Gaskill’s  cylinder 
function: 

=  Ziyi 

Thus  fv{v)  is  an  average  function  over  the  entire  entrance  aperture  of  each 
mini-reflector.  But  because  the  retro-material  contains  many  thousands  of  mini¬ 
reflectors,  rather  than  compute  individual  lines  of  sight  through  the  retro-material 
geometry,  it  is  more  reasonable  to  apply  /y('u)  as  a  statistical  quantity  for  any 
radiation  reaching  the  retro-material  plane  b\.z  =  Z.  In  effect,  then,  fv{v)  acts  to 
reduce  the  reflected  radiance  depending  on  the  angle  of  incidence  of  the  received 
radiance  relative  to  a  normal  vector  to  the  plane  of  the  retro-material.  And  because 
fv{u)  has  a  multiplying  effect  on  the  radiance  magnitude,  upon  Fourier  transfor¬ 
mation,  the  result  is  a  convolution  operation  in  angular  frequency  that  is  applied 
to  the  return  beam.  This  angular  filtering  also  tends  to  reduce  the  importance  of 
large  angle  scattering  because  energy  arriving  from  large  angles  is  reduce  greatly 
in  magnitude  on  retroreflection. 


In  addition  to  this  magnitude  modulation  as  a  function  of  incidence  angle,  a  spatial 
truncation  occurs  due  to  the  finite  size  of  the  retro-material.  We  assume  the  retro- 
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material  to  be  circular  in  shape,  with  a  diameter  Dr,  and  oriented  perpendicular 
to  the  2  axis.  Let  the  center  of  this  reflective  material  be  located  off  axis  at  point 
f R  in  ihez  =  Z  plane.  Due  to  scenario  considerations  we  will  place  this  off  axis 
position  somewhere  along  the  x  axis.  Mathematically,  this  reflector  region  could 
be  represented  by  another  filter  function, 

/w(f)  =  cyl 

However,  from  a  practical  consideration,  the  size  of  the  reflector  region  will  likely 
be  small.  Thus,  relative  to  the  overall  spread  in  the  beam  itself  and  the  distance 
off  axis,  we  will  consider  the  directional  and  positional  properties  of  the  incident 
beam  to  be  approximately  constant  over  the  surface  of  the  retroreflector.  We  can 
then  just  as  easily  compute  the  reflected  radiation  via  use  of  a  5  function  and  an 
equivalent  reflector  area. 

Finally,  diffraction  effects  at  the  retroreflector  will  also  influence  the  angular  dis¬ 
tribution  of  the  retroreflected  energy.  While  the  total  retroreflector  area  has  a  di¬ 
ameter  Dr,  the  retro-material  is  actually  composed  of  thousands  of  mini-reflectors 
with  effective  aperture  radii  of  only  as  in  figure  12.  We  can  model  the  diffrac¬ 
tion  spread  of  the  return  beam  due  to  the  area  of  these  mini-reflectors  as  an  angular 
point  spread  function.  However,  crude  experiments  indicate  the  actual  beamspread 
may  be  1.5-3  times  larger  than  the  diffraction  limited  results.  In  frequency  space 
this  point  spread  function  results  in  an  additional  modulation  transfer  function 
applied  to  the  return  beam. 

These  influences  of  the  retroreflective  surface  in  spreading  the  retroreflected  beam 
are  advantageous  since  the  system  envisioned  would  locate  the  observation  cam¬ 
era  at  a  distance  offaxis  from  the  laser  provideing  the  source  illumination.  Some 
beam  spread  thus  increases  the  amount  of  energy  that  can  be  received  by  the  off- 
axis  camera.  Insufficient  spreading  would  result  in  no  received  energy,  while  too 
much  spread  would  likely  result  in  further  blurring  of  the  received  image. 

Combining  these  considerations,  it  may  be  possible  to  approximately  represent 
the  return  beam  as  a  Gaussian  beam,  with  the  return  path  adjusted  in  such  a  way 
that  the  retroreflector  appears  on  axis  and  the  reflected  “Gaussian  beam  appears 
to  be  travelling  along  the  x  =  y  =  0  axis.  However,  the  receiver  camera  would 
need  to  be  offset  from  this  new  optical  axis.  To  handle  these  considerations  we 
would  need  to  assess  the  direction  associated  with  the  main  retroreflecting  beam 
centroid,  call  it  ur.  Let  fx  be  the  offset  position  of  the  receiver  at  the  2;  =  0  plane 
from  the  source  laser  at  fi  =  0.  Let  the  main  direction  of  incident  energy  arriving 
at  the  retroreflector  position  be  ujr.  Then  the  direction  of  the  main  radiance  reflect¬ 
ing  back  from  the  retroreflector  along  the  return  path  would  be  travelling  in  the 
direction  associated  with  the  triplet  {—ujr,  — !)•  Projecting  this  main  propagation 
centroid  back  to  the  z  =  0  plane  we  would  find  an  intercept  at  the  point 

fc  =  ^R  —  ZujR.  (Ill) 
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Let  us  then  define  a  new  coordinate  system  in  which  we  relocate  fc  onto  the  main 
optical  axis  x'  —  y'  =  0.  In  this  new  system  the  retroreflected  energy  appears  to 
be  travelling  similar  to  the  original  beam.  The  only  changes  that  are  needed  are 
to  assess  the  characteristics  of  the  new  beam  and  to  displace  the  position  of  the 
receiver: 

f^x  =  rx-fc.  (112) 

Having  made  this  change  in  geometry  all  that  remains  in  describing  the  return 
radiation  field  is  to  examine  the  width  of  the  Gaussian  beam.  We  then  consider 
the  propagation  problem  along  the  return  path  in  chapter  4. 

To  describe  the  additional  spreading  of  the  return  beam  due  to  the  retroreflector, 
we  begin  by  assuming  we  can  characterize  the  incident  beam’s  angular  stmcture 
as  a  Gaussian  about  some  mean  direction  wr.  We  similarly  assume  we  can  charac¬ 
terize  the  nondiffraction-limited  effects  of  the  retroreflector  aperture  as  a  Gaussian 
as  well.  These  effects  are  then  combined  to  generate  a  third  Gaussian  form  for  the 
output  from  the  retroreflective  material  about  the  mean  retro-direction  which  we 
have  remapped  to  ujr  =  0.  Let  vi  represent  the  angular  width  of  the  incident 
radiation,  and  let  V2  represent  the  spreading  effect  of  the  aperture  on  a  monodirec- 
tional  incident  beam.  Then  we  characterize  the  total  output  beam  angular  width 
parameter,  v,  by 

v‘^  =  vf  +  vl  (113) 

In  evaluating  the  reliability  of  this  model  we  note  that  when  the  fog  is  thick,  there 
will  be  significant  effects  on  the  spreading  of  the  beam  with  location  and  propa¬ 
gation  direction  on  the  return  passage  through  the  fog.  It  is  therefore  not  critical 
to  exactly  model  the  nature  of  the  initial  angular  structure  of  the  return  beam.  Of 
much  more  critical  interest  is  that  we  have  captured  the  positional  offset  of  the 
retroreflector  from  the  optical  axis  by  the  use  of  vr. 

The  remaining  task  is  to  evaluate  V2.  An  approximation  of  this  quantity  can  be 
made  based  on  the  observation  that  there  are  approximately  60  mini-reflectors 
per  inch.  We  thus  have  Xm  «  210  //m.  For  an  aperture  of  diameter  D  the  far 
field  diffraction  pattern  will  have  a  zero  at  an  angle  of  l.22\/D.  The  amplitude  of 
this  pattern  has  a  form  called  an  Airy  disc,  or  in  Gaskill’s  terminology,  a  sombrero 
function.  In  either  case,  because  the  light  passing  through  the  mini-reflectors  actu¬ 
ally  passes  through  two  circular  apertures,  it  is  believed  that  the  actual  pattern  for 
the  amplitude  follows  the  square  of  a  sombrero  function.  Therefore,  the  irradiance 
pattern  generated  by  each  mini-reflector  is  modeled  by  the  sombrero  function  to 
the  fourth  power.  This  pattern  may  be  approximated  by  a  Gaussian  function  with 
an  angular  width  parameter  of  approximately  V2  =  OSlX/vm,  where  V2  appears 
as  a  width  parameter  in  the  irradiance  form:  X(r',  z')  oc  exp[— 2r'^/(y^  v^)].  To 
use  this  result  we  will  also  need  to  evaluate  the  angular  properties  of  the  incident 
energy  to  determine  a  central  peak  in  the  propagation  direction  and  to  determine 
the  angular  deviation  that  results  in  a  reduction  of  exp(-2)  from  the  central  peak 
value.  The  angular  value  of  this  result  then  becomes  vi  and  can  be  used  to  compute 
V  in  equation  111. 
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4.  Return  Path  Propagation 


Figure  14  considers  the  overall  geometry  of  the  problem  we  are  modeling.  Though 
we  appear  to  have  placed  the  retroreflector  along  the  optical  axis  of  the  source, 
this  is  not  necessary.  The  key  point  is  that  the  total  path  need  not  be  considered 
folded.  Rather,  the  path  can  be  viewed  as  three  separate  segments  of  one  continu¬ 
ous  path.  The  first  portion  leads  from  the  laser  source  to  the  entrance  pupil  of  the 
retroreflector.  The  second  stage  (treated  in  chapter  3)  passes  the  beam  through  the 
retroreflector.  The  third  segment  includes  the  return  path,  which  is  a  mirror  image 
of  the  outbound  propagation  path. 


Retroreflector 


Laser  Source 

h 


Fog  Layer 

I  r  I  m  n  m m  *  ■  I ! !  *  H 

2m  — J— — Z  -  Zm 


.  Wing  Camera 


- - 

- Z-7.^ 

:::  1^ 

’l—Zm-A 

- ^Forward  Path - -  - Return  Path - - 

Figure  14.  Geometry  of  the  folded  path  problem  where  we  treat  the  return  path  as  a  reversed  version  of 
the  forward  propagation  path. 


In  studying  the  return  portion  of  the  propagation  path  we  note  that  the  solu¬ 
tion  found  for  the  outbound  propagation  path  was  radially  symmetric.  The  return 
beam  will  also  share  this  symmetry  property.  In  our  problem,  the  laser  source  is 
mounted  somewhere  on  the  fuselage  or  wing  of  an  aircraft.  Receiver  cameras  will 
be  mounted  on  each  of  the  wings.  Hence,  if  the  laser  and  two  cameras  are  mounted 
along  an  axis  of  the  wing,  and  the  aircraft  is  flying  level  to  the  terrain  generally  in 
the  direction  toward  the  landing  field,  it  is  likely  that  at  least  to  a  first  approxima¬ 
tion  the  cameras,  laser,  and  retroreflector  can  all  be  considered  to  be  lying  within 
a  single  plane.  And  somewhere  along  the  line  containing  the  laser  source  and  the 
cameras  the  return  signal  will  be  a  maximum.  This  will  simplify  the  subsequent 
analysis  since  we  will  only  need  to  look  along  a  single  line  (call  it  the  x  axis)  in 
determining  the  position  of  the  maximum  return  signal. 

Based  on  the  discussion  in  chapter  3,  we  will  assume  the  receiver  position  has 
already  been  adjusted  to  the  location  f^.  It  should  also  be  understood  that  in  this 
chapter  we  are  going  to  be  dealing  with  the  return  path.  Nevertheless,  to  avoid 
complications  with  notation  we  are  not  going  to  adopt  a  completely  different  sys¬ 
tem  to  denote  the  different  axes  of  propagation  used  in  analyzing  the  return  prop¬ 
agation.  It  will  have  to  be  assumed  that  when  using  such  terms  as  x,  y,  z,  f,  and  u. 
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in  this  chapter  that  we  are  nevertheless  talking  about  a  different  problem  than  was 
considered  in  either  chapters  2  or  3.  We  will  see  that  these  problems  are  still  here 
to  haunt  us  in  chapter  5  where  we  treat  radiation  propagating  in  both  hemispheres 
at  once,  but  we  will  leave  that  as  a  separate  problem.  For  this  chapter  we  simply 
adopt  the  same  terminology  as  was  used  in  chapter  2,  but  apply  it  to  the  return 
path. 

We  also  mentioned  at  the  close  of  chapter  3  that  we  neededto  apply  a  factor  v, 
which  represented  the  total  effects  of  scattering  on  the  outbound  path  and  diffrac¬ 
tion  interaction  with  the  retro-material.  In  chapter  2  we  characterized  the  original 
Gaussian  beam  emitted  by  the  source  laser  expressed  in  radiance  form  in  equa¬ 
tion  (72).  For  the  return  path  we  parameterize  the  return  beam  by, 

(114) 


I^{f,  u),  Z)  = 


2Pr 

2  Z2 


-nv 


exp 


-2r^ 


where  Pr  is  the  total  reflected  power  of  the  return  beam.  Note  that  if  Vi  was  the 
effective  angular  spread  of  the  beam  arriving  at  the  retro-material,  and  Rr  was 
the  radiance  arriving  at  the  material  in  the  peak  direction  of  incidence,  then  the 
arriving  radiance  as  a  function  of  angle  would  be 


Z)  w  Rr  exp 
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where  we  assume  the  incident  irradiance  is  roughly  constant  over  the  surface  of 
the  retroreflector.  Integrating  with  angle  and  multiplying  by  the  surface  area  of  the 
retroreflector,  we  have 

(116) 

where  Ar  represents  the  effective  area  of  the  retroreflector.  This  factor  may  take 
into  account  the  foreshortening  of  the  reflective  surface  if  it  is  not  perpendicular  to 
the  incident  beam  and  the  effects  of  off-axis  angle  of  arrival  of  the  incident  beam 
to  the  retroreflective  material. 


Once  the  initial  return  beam  is  postulated,  we  may  proceed  to  evaluate  the  net 
returned  beam  passing  through  the  scattering  media  similar  to  the  method  used  in 
chapter  2.  The  key  difference  here  is  that  the  fog  aerosol  will  occupy  the  initial 
portion  of  the  path  whereas  previously  it  occupied  the  final  segment  of  the  path. 

In  solving  for  the  return  beam  we  first  note  that  due  to  our  change  of  geometry 
we  have  the  same  basic  problem  to  solve  for  both  the  vacuum  and  unscattered 
problems.  For  completeness,  we  include  these  results  here.  First,  for  the  vacuum 
solution,  given  the  slight  changes  involved  with  the  introduction  of  v,  we  have 


/^(k,  fj,  Z)  =  Pr  exp 


^2^2  ^2  p^2^z) 
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We  then  introduce  a  simplify  A  form, 
A(z)  ~  K^-rj/Z 


{fjz  +  Zk) 


(118) 
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such  that  we  can  express  this  result  as 


/^(k,  ff,  Z)  =  Pr  exp 


TT^  T]q 
2 


(119) 


where  zero  range  is  now  viewed  as  lying  at  the  retroreflector  plane. 

The  unscattered  solution  then  follows  simply  as 

exp  [-t(Z)]  ,  t{Z)  =  aQZ.  (120) 


It  then  remains  to  evaluate  the  scattered  solution.  This  form  is  slightly  different 
from  the  forward  propagation  problem  because  Sr  depends  on  aerosol  that  occu¬ 
pies  the  first  portion  of  Ae  optical  path: 


Zf  ^ 

Sr  =  exp  ^  w  a  J  Pf  {{Vz  +  k[Z  -  z]\)  dz 
0 


=  n  ^ 


^  {.Qm  “b  (1  Q')Qm  ~b  I^Pn 


.  (121) 


However,  as  expected,  the  resulting  equation  for  Sr  generates  the  same  results 
when  the  aerosol  occupies  the  full  path  (Q  =  1). 

The  same  techniques  are  used  to  evaluate  the  inverse  transform  of  this  result  as 
were  discussed  in  chapter  2  for  handling  the  forward  propagation  to  the  retrore¬ 
flector  plane. 


5.  Diffuse  Radiation 


In  previous  chapters  we  considered  the  forward-scattered  radiation  propagating 
both  from  the  laser  source  to  the  retroreflector  and  from  the  retroreflector  back 
to  one  of  the  wing  cameras.  In  addition  radiation  may  scatter  from  the  fog  vol¬ 
ume  itself  and  be  detected  by  one  of  the  wing  cameras.  This  scattered  radiation 
appears  as  a  noise  source  that  must  eventually  be  distinguished  from  the  return  en¬ 
ergy  propagating  back  from  the  retroreflector.  Several  strategies  were  mentioned 
for  removing  this  additional  radiation.  The  simplest  solution  involves  image  pro¬ 
cessing  to  discriminate  the  high  angular  resolution  return  signal  from  the  diffuse 
backscattered  radiance  arising  from  scattering  off  the  fog  alone. 

To  determine  the  difficulty  of  the  discrimination  problem,  the  diffuse  scattering  ra¬ 
diance  field  must  be  modeled.  The  initial  Gaussian  beam  forms  the  energy  source 
for  this  radiation  field.  However,  unlike  the  forward  scattering  problem  we  need 
significantly  more  information  to  perform  this  modeling  task.  To  be  able  to  model 
the  backscattered  radiance  field  we  will  need  to  know  details  of  the  forward  scat¬ 
tered  field  at  all  locations.  Determining  this  information  using  the  full  resolution 
phase  function  does  not  appear  practical.  Hence,  we  use  a  simplified  version  of 
the  forward  scattered  field.  Second,  we  will  ignore  the  small  amount  of  energy 
retroreflected  from  the  comer  reflector.  Relative  to  the  energy  scattering  within 
the  fog  cloud  this  amount  will  be  minimal. 

5.1  Diffuse  Radiation  Scale  Transform 

To  perform  this  analysis,  the  total  diffuse  radiation  field  will  be  divided  into  two 
portions:  (1)  Idf  flowing  in  the  positive  z  direction  >  0);  (2)  hb  flowing  in 
the  negative  2  direction  (f22  >  0).  Previously  we  used  Q  to  represent  the  x  and  y 
components  of  the  propagation  vector  in  the  positive  z  direction  only.  Here,  we 
will  also  use  u  to  denote  propagation  in  the  negative  z  direction  but  will  use  the 
following  convention: 


~  — tt’y,  1)  •  (122) 

Hence,  Idf  and  Idb  radiances  assigned  equal  u  arguments  represent  energy  prop¬ 
agation  in  opposite  directions.  This  property  will  be  useful  when  evaluating  the 
scattering  integral  in  the  RTE. 

In  addition  to  Idf  flowing  in  the  positive  z  direction,  we  will  introduce  the  quantity 
If{f,  u),  z)  representing  energy  that  is  propagating  in  nearly  the  forward  direction. 
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Ijf  and  Idb  thus  represent  generalized  diffuse  fields  with  low  angular  resolution, 
while  the  //  field  constitutes  both  unscattered  energy  from  the  laser  source  and 
forward-scattered  energy  in  a  small  solid  angular  region  about  the  forward  direc¬ 
tion.  To  model  If  we  use  a  considerably  simpler  RTE  than  was  used  in  chapter 

where  the  quantity  Af  represents  the  probability  of  scatter  in  the  “forward”  di¬ 
rection.  This  quantity  can  be  viewed  as  one  of  the  Am  weighting  factors  used 
in  chapter  2.  The  remainder  of  the  phase  function  can  then  be  modeled  using  a 
shorter  Gaussian  expansion.  In  particular,  we  will  use  the  same  expansion  terms 
as  are  contained  in  table  1,  except  that  the  first  two  terms  have  been  removed.  The 
corresponding  A/  is  thus  the  sum  of  the  first  two  Am  terms  in  table  1  and  has  the 
value  0.4567.  Because  oi  and  0-2  are  significantly  greater  than  the  terms,  we  are 
justified  in  assuming  a  zero  width  for  the  forward  scattering  portion. 

This  form  for  the  radiative  transfer  of  the  nondiffusely  scattered  radiation  assumes 
the  forward  scattered  energy  is  scattered  into  a  cone  so  small  that  even  scattered 
energy  can  be  assumed  to  be  propagating  virtually  along  the  path  of  the  original 
unscattered  energy.  To  solve  the  If  equation  we  can  perform  a  scale  transforma¬ 
tion  on  the  extinction  coefficient  (McKellar  and  Box,  1981).  This  approach  in¬ 
volves  the  replacement  of  the  original  extinction  coefficient  by  the  reduced  quan¬ 
tity: 

0  —  —  AfZj)  a.  (124) 

Through  this  replacement  the  RTE  for  If  is  reduced  to  the  unscattered  RTE  with 
cr  as  the  extinction  coefficient. 

d5.^  +  ^  +  5//=0,  (125) 

or  oz 

which  has  the  solution. 

If  =  I^  exp  [-t(2:)]  ,  (126) 

where  I^  is  the  vacuum  solution  introduced  in  equation  (72),  and 


t{z)  =  J  ads.  (127) 

0 


A  scaling  transformation 
which 


also  entails  a  change  to  the  single  scattering  albedo  in 


{1-Af)w 
(1  —  Afw)  ’ 


(128) 


such  that  the  product  zaa  =  (1  -  Af)'cja  appears  as  a  scattering  coefficient  (c.f. 
Tofsted  and  O’Brien,  1998). 


The  next  task  is  to  write  RTE’s  for  both  the  forward  and  backward  oriented  diffuse 
radiance  fields.  To  write  these  equation,  however,  we  must  recognize  that  scatter¬ 
ing  occurs  between  radiances  flowing  into  opposite  hemispheres.  But  since  the 


phase  functions  are  generally  written  in  terms  of  scattering  angles  measured  by 
Au,  we  must  make  allowances  for  scattering  into  opposite  hemispheres  through 
the  use  of  two  different  phase  functions.  The  first, 

M 

PdfiAoj)  =^A'^al,  exp  (-a^  Aa;^)  ,  (129) 

m=3 

represents  scattering  in  the  forward  direction  (within  the  same  hemisphere  into 
which  the  scattering  energy  was  originally  travelling).  The  second, 

Pdb(  Aw)  =  exp  {-Pi  Aw^)  ,  (130) 

n=l 

represents  the  scattering  of  energy  travelling  in  the  direction  of  one  hemisphere 
into  the  opposite  hemisphere.  It  both  these  cases  the  A'^  and  quantities  have 
been  adjusted  from  the  values  listed  in  table  1  by  the  following  adjustments: 

=  Am!  {I  -  Af)  and  B'^  =  5„/(l  -  Af).  The  reason  for  these  adjustments  is 
that  due  to  the  scale  transformation,  the  diffuse  scattering  terms  represent  100  per¬ 
cent  of  the  scattering  in  the  scale  transformed  equations.  But  the  sum  of  the  scat¬ 
tering  components  (E^m  +  E-^n)  now  only  accounts  for  (1  -  A/)  probability 
of  scattering.  To  renormalize  the  phase  functions  we  need  to  divide  by  this  factor. 


5.2  Dual  Direction  Radiance  Equations 

We  can  then  write  the  coupled  differential  equations  describing  the  propagating 
forward  and  backward  diffuse  radiance  fields  as 

w  •  ^  ^  +  aldf  =  wa  jjPdf{Au)  (/}  -I-  du' 

—  OO 
cx» 

+  w  a  jJPdb{Auj)  (131) 

—  OO 

oo 

-w  •  +  0^1  =  wa  jjPdb{Auj)  [I'j  -k  duj', 

—  OO 
OO 

-k  ^Pd/(Aw)4dw',  (132) 

—  OO 

where  Aw  =  |w  -  w'|  as  usual.  The  terms  /},  and  refer  to  radiances 
propagating  in  directions  associated  with  the  vector  w'.  Note  that  the  scattering  of 
radiance  out  of  the  opposite  hemisphere  is  centered  on  the  backscattering  direc¬ 
tion.  Also  note  that  the  directional  derivative  signs  in  the  Idb  equation  are  switched 
to  represent  the  opposite  direction  of  propagation  for  the  backward  propagating 
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radiance  field.  Evaluating  both  these  diffuse  radiance  fields,  and  coupling  the  re¬ 
sults  of  this  analysis  to  the  results  of  computations  based  on  analysis  contained 
in  chapters  2-4,  we  create  a  complete  picture  of  the  radiation  flowing  through  the 
fog  cloud. 

To  solve  this  pair  of  equations,  we  first  proceed  as  previously  by  taking  the  fourth 
order  Fourier  transforms  (in  r  and  u)  of  both  equations: 

-/c  •  ^  -f  ^  -I-  o-  idf  =  w5-  ^Pdf  j^/d/  +  //  +  Pdb  ,  (133) 

-t-K  •  ^  ^  -h  a  idh  =  '^Sr  jp*  \^df  +  //  +  Pdf  /dftj  •  (134) 

The  physical  geometry  of  the  problem  to  be  solved  can  thus  be  stated  as  follows: 
We  assume  that  the  leading  edge  of  the  cloud  begins  at  some  distance  from  the 
laser  source.  This  distance  was  called  Zf  in  chapters  2  and  4.  Let  us  then  assume 
that  the  cloud  ends  at  distance  Zb.  The  length  of  the  cloud  is  then  Zc  and  at 
distance  Z  embedded  within  the  cloud  is  the  retroreflector. 

Since  no  forward  scattering  may  occur  before  the  forward  edge  of  the  cloud  we 
must  have 

idf{Zf)  =  0.  (135) 

Likewise,  there  must  be  no  backward  scattering  beyond  the  farthest  edge  of  the 
cloud: 

idbiZb)  =  0.  (136) 

We  may  also  attempt  to  simplify  the  system  of  equations  by  adding  and  subtracting 
eqs.  (131)  and  (132)  from  one  another.  Combining  terms  we  obtain  the  equation 

/l -t- 5-/2  =  rocrPu  12  + if  ,  (137) 

/2  +  ail  =  wa  Ps  ii  +  if  ,  (138) 

where  ....  ... 

il  =  idf  +  idbi  ^2  =  idf  ~  idbi  (139) 

Ps  =  Pdf  +  Pdb,  Pd  =  Pdf  —  Pdb-  (140) 

For  this  system  we  have  boundary  conditions, 

il  =  — /2>  Z  =  Zf',  il  =  +/2)  ^  —  ^b-  (141) 

Given  these  boundary  conditions  we  consider  the  general  solution  to  this  set 
of  equations.  The  first  question  we  need  to  ask  is  whether  the  implicit  solution 
method  introduced  in  chapter  2  can  be  used  in  solving  this  problem.  It  would 
seem  that  the  implicit  method  should  work  individually  on  each  pair  of  equations. 


'  ^  d  d' 
^  dr]^  dz 

^  dfj  dz 
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Hence,  we  have  two  analyses  similar  to  the  development  leading  to  eqs.  (63),  (64), 
and  (70).  In  particular,  we  generate  the  same  relationship  in  both  eqs.  (135)  and 
(136)  to  describe  the  behavior  of  t?  as  a  function  of  k  and  z: 


J7(z)  =  ffo-KZ. 


(142) 


The  remaining  relationships  provide  differential  equations  for  the  radiances  in 
terms  of  distance; 


dz 


—a 


J2  - 


(143) 


dz 


—a 


(144) 


This  system  of  equations  can  be  solved  through  use  of  the  Newton-Raphson  method. 

Solution  involves  determining  a  starting  value  for  1 1  at  Zf  that  when  propagated 
through  to  the  far  end  of  the  cloud  at  range  Zb  produces  the  correct  boundary  con¬ 
dition  listed  in  equation  (139).  Since  there  is  no  change  in  either  variable  when  d 

is  zero,  and  d  is  zero  prior  to  the  beginning  of  the  cloud,  the  values  of  both  1 1  and 
/2  at  zero  range  are  also  their  values  at  the  forward  edge  of  the  cloud. 


The  Newton-Raphson  method  can  be  used  by  setting  the  function  we  wish  to  find 

a  zero  for  as  7i  -  /2  evaluated  at  Zb  as  a  function  of  /i  at  Z/ .  The  means  of  evaluat¬ 
ing  the  two  output  parameters  at  the  far  end  of  the  cloud  involves  implementation 
of  a  numerical  integration  procedure.  Press  et  al.  (1992)  provide  listings  of  proce¬ 
dures  to  accomplish  these  tasks,  including  routines  rtsaf  e  to  do  a  combination 
Newton-Raphson  and  bisection  procedure,  ode  in  t  to  integrate  a  set  of  ordinary 
differential  equations,  and  codes  rkqs  and  rkck  to  do  quality-controlled  fourth- 
order  Cash-Karp  Runge-Kutta  integration  of  the  differential  equations. 

Once  the  Newton-Raphson  approach  has  been  used  to  evaluate  Ji,  we  can  evaluate 
Idb  at  2  =  0  using  the  same  inverse  transform  techniques  described  in  sections 
2.7.5  and  2.8.  The  only  difference  is  in  the  interpretation  of  the  outputs:  While 
the  f  system  will  be  the  same,  the  to  sense  of  the  incoming  radiance  field  will  be 
the  opposite  of  the  sense  of  the  outbound  radiances  of  the  same  magnitude  and 
direction. 


We  recover  the  original  diffuse  fields  through  the  equations 


id! 


(145) 


With  this  result  we  have  shown  the  means  of  computing  all  the  fields  and  quanti¬ 
ties  of  interest.  We  leave  the  analysis  of  these  results  to  a  future  publication. 
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6.  Conclusions 


This  document  develops  a  theoretical  basis  for  the  problem  of  dual  path  propaga¬ 
tion  through  scattering  aerosol  advection  fog.  To  solve  this  problem  our  approach 
has  been  rigorous  —  we  have  considered  both  propagation  through  a  fog  layer  to 
a  retroreflector  and  return  propagation  through  the  fog.  We  have  also  considered 
the  issue  of  scattering  within  the  fog  layer  in  both  forward  and  backward  direc¬ 
tions.  A  recasting  of  the  small  angle  approximation  has  also  been  accomplished 
as  well  as  a  complete  elucidation  of  the  methods  employed  by  Smirnov  to  solve 
partial  differential  equations,  and  this  approach  has  been  applied  to  the  propaga¬ 
tion  problem  at  hand.  In  addition,  we  have  considered  means  useful  in  reducing 
the  inverse  transform  problem  from  a  fourth  order  inverse  Fourier  transform  pro¬ 
cess  to  a  third  order  calculation  involving  Bessel  functions.  This  method  greatly 
simplifies  the  inverse  transform  process  since  it  guarantees  that  one  will  generate 
a  real-valued  transform  function  as  a  result.  We  then  considered  the  process  of 
retroreflection  itself  and  provided  a  model  useful  in  approximating  the  effects  of 
retroreflection  from  two  commercially  available  retroreflective  materials.  Lastly, 
we  considered  the  calculation  of  diffuse  radiation  scattering  within  the  fog  cloud 
volume.  While  complicated,  these  procedures  provide  an  accurate  means  of  evalu¬ 
ating  the  effects  of  propagation  within  fog  aerosols  to  assess  a  navigation  through 
fog  system. 
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